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Section  1 
INTRODUCTION 


The  statistical  approach  to  the  solution  of  radiant  energy  transfer  among 
surfaces  of  a  system,  commonly  known  as  the  Monte  Carlo  Method,  makes 
it  possible  to  attack  difficult  problems  which  are  not  amenable  to  conven¬ 
tional  analytical  or  computer  methods  such  as  Fourier  Transform  tech¬ 
niques.  Monte  Carlo  methods  represent  radiant  energy  flux  statistically 
as  rays  which  can  be  traced  through  the  system  according  to  a  set  of 
probability  distribution  functions.  Although  the  ray  representation  of  the 
radiant  energy  is  similar  to  that  of  a  photon  character;_ation,  it  is  treated 
in  this  method  as  merely  a  statistical  representation. 

A  distinct  disadvantage  of  the  method  is  the  extremely  large  number  of 
calculations  necessary  to  simulate  a  fairly  complicated  system  with  a 
reasonable  degree  of  accuracy.  With  the  advent  of  high  speed  digital  com¬ 
puters,  however,  it  has  become  a  feasible  tool  in  a  wide  variety  of  problems. 

With  the  general  radiation  analysis  in  mind,  the  GUERAP  (General  Unwanted 
Energy  Rejection  Analysis  Program)  was  written  specifically  for  tracing 
radiant  energy  through  systems  with  extremely  high  attenuation.  Since  the 
number  of  rays  necessary  to  analyze  a  system  normally  is  directly  pro¬ 
portional  to  its  attenuation,  the  "brute  force"  Monte  Carlo  method  becomes 
impractical  in  tracing  rays  through  systems  with  high  attenuations,  e.  g. , 
>106. 

Many  rays  are  required  for  a  high  attenuation  system  because  all  but  a 
handful  of  the  rays  generated  are  lost  in  the  system,  without  contributing 
anything  to  the  net  energy  at  points  of  interest,  e.g. ,  a  detector.  The 
situation  can  be  drastically  improved  if  rays  are  directed  to  follow  pre¬ 
determined  paths  toward  the  points  of  interest. 

This  is  accomplished  by  dividing  the  angular  probability  density  distribu¬ 
tion  function  of  an  emerging  ray  into  multiple  sections.  Each  section  is 
represented  by  a  ray  whose  direction  is  randomly  selected  within  the 
angular  confine  of  the  section.  The  total  energy  of  the  original  ray  is 
distributed  among  the  several  rays  according  to  their  "share"  of  the  total 
probability  distribution  function.  This  method  is  commonly  known  as  the 
Importance  Sampling  or  Expected  Value  Technique  (Reference  1). 
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Diffraction  of  radiant  energy  as  it  passes  through  an  aperture  is  generally 
regarded  as  a  manifestation  of  the  wave-like  nature  of  radiation.  The 
theory  then  dictates  that  the  diffraction  pattern  behind  the  aperture  can 
only  be  obtained  by  simultaneously  considering  the  total  wave  in  the  sys¬ 
tem,  thus  making  it  impossible  to  rigorously  include  diffraction  effects  by 
the  Monte  Carlo  ray  tracing  approach. 

A  semi-empirical  diffraction  model  is  adopted  in  this  program  that  ansunu  » 
an  angular  probability  density  function  of  a  ray  as  it  passes  through  a  small 
aperture  (Reference  2).  The  model  makes  use  of  the  prediction  by  the 
Uncertainty  Principle  (Reference  4)  which  allows  the  ray  to  be  b*nt  as  a 
result  of  the  interaction  between  the  ray  and  the  aperture  edge.  The  model 
is  successfully  tested  by  comparing  results  of  several  simple  systems 
with  classical  solutions. 

To  analyze  radiant  energy  transfer  in  a  system,  the  user  first  represents 
the  system  with  a  set  of  surfaces,  expressible  in  terms  of  simple  equa- 
i  ons.  A  series  of  basic  surfaces  that  can  be  used  for  such  purpose  arc 
presented  in  Section  2.  Although  it  is  sometimes  necessary  to  supplement 
by  inputting  surfaces  in  terms  of  general  equations,  these  basic  surfaces 
are  sufficient  for  most  systems  commonly  encountered.  We  then  proceed 
to  present  the  procedure  of  simulating  the  systems. 

In  Section  3,  details  are  given  for  various  radiation  models.  These  include 
emission,  absorption,  reflection,  transmission,  refraction  and  diffraction. 
The  total  energy  of  a  ray  interacting  with  a  surface  is  represented  by  an 
angular  distribution  function  ior  each  model. 

The  importance  sampling  technique  is  briefly  described  in  Section  4.  For 
the  usage  of  the  technique,  the  user  is  referred  to  the  procedure  of  assign¬ 
ing  important  surfaces  given  in  Section  2  as  well  as  Section  4. 

For  users  who  wish  to  gain  better  understanding  of  the  program.  Section  5 
presents  a  brief  description  of  the  program  structure  and  the  purpose  of 
each  subroutine.  Flow  charts  are  given  for  the  main  program  and  the 
commanding  subroutines. 

The  complete  set  of  input  parameters  is  listed  in  Section  6.  Each  param¬ 
eter  is  followed  by  a  brief  definition  and  the  page  number  for  reference. 

An  experienced  user  of  the  program  can  complete  a  set  of  input  data  with 
the  help  of  this  section  alone. 


1  -? 


3366-008 


Most  of  the  output  listing  and  the  error  messages  are  self-explanatory. 
However,  the  user  can  refer  to  Sections  7  and  8  for  more  detailed  explana¬ 
tions  lacking  in  the  output  listing. 

Finally,  a  series  of  examples  are  presented  in  Section  9.  Experiences  on 
the  part  of  this  writer  indicate  that,  in  learning  the  usage  of  a  new  program, 
it  is  sometimes  easier  to  follow  the  lead  of  examples.  For  this  reason, 
the  examples  are  presented  for  problems  with  different  degrees  of  complex¬ 
ity.  As  many  features  as  deemed  necessary  to  demonstrate  the  common 
usage  of  the  program  are  included  in  the  examples. 
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Section  2 

SYSTEM  CONFIGURATION 


The  basis  of  the  Monte  Carlo  method  is  in  tracing  the  probable  paths  of  an 
energy  bundle  as  it  interacts  with  the  optical  system.  The  present  program 
is  concerned  primarily  with  radiant  energy  transfer  through  non-participating 
media.  Thus,  interactions  occur  only  when  a  ray  is  intercepted  by  an  optical 
surface  of  the  system. 

In  order  to  detect  such  inter  ceptions,  it  is  necessary  to  represent  the  sur¬ 
faces  by  some  mathematical  model,  preferably  in  the  form  of  analytical 
equations  that  define  the  surfaces  and  their  boundaries.  A  series  of  simple 
surfaces  have  been  developed  as  the  basis  of  a  general  optical  system. 

While  these  surfaces  are  sufficient  for  the  simulation  of  the  systems  en¬ 
countered  during  the  development  of  the  program,  some  additional  basic 
surfaces  will  undoubtedly  be  added  as  needs  arise. 

The  basic  surfaces  with  their  standard  boundaries  (the  boundaries  defined 
by  the  basic  parameters)  will  be  first  described.  The  boundaries,  which 
we  refer  to  as  constraints,  are  represented  by  a  set  of  constraint  equations. 
Frequently,  a  surface  requires  a  set  of  constraints  other  than  the  standard 
set.  Thus,  attention  will  be  given  to  the  method  of  prescribing  these 
special  constraints. 

When  a  surface  cannot  be  represented  by  any  of  the  basic  surfaces,  the 
coefficients  of  the  surface  equation: 

C,x2  ♦  C2Y2  ♦  C3Z2  ♦  C4XY  *  CgXZ+  CgYZ  ♦  C,X  ♦  CgY  *  CgZ  ■  C,0 

(2-1) 

can  be  read  in  through  the  array 

COEF(J.I)  =  Cj  J  ■  1,  2 . 10 

where  I  is  the  surface  number.  The  degree  of  the  equation  is  specified  by 
the  array  IDGREE  in  input  data  set  NAME2.  Namely,  IDGREE(T)  is  set  at 
1  if  the  first  six  coefficients  in  the  surface  equation  are  zero.  Otherwise, 
it  is  set  equal  to  2.  The  set  of  constraints  for  the  surface  is  specified  as 
special  constraints  (see  Subsection  2.2). 
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2.  1  BASIC  CONFIGURATION 

2.1.1  Plane 

A  plane  bounded  by  two  concentric  circles  as  shown  in  Figure  2-1  is  used 
as  the  basic  plane  configuration.  It  is  described  by  the  coordinates  of  the 
center  point  (Xc,  Y( ,  Zc),  the  unit  normal  vector  (a,  (J,  y)  and  the  radii  of 
constraints  Rj  and  R2.  When  it  is  possible,  the  unit  normal  should  point 
into  the  region.  The  parameters  in  this  sequence  are  stored  as  PTLCTR 
of  the  surface. 

In  subroutine  PLANE  coefficients  of  the  surface  and  the  constraint  equations 
in  terms  of  the  system  coordinates  X-Y-Z*  are  calculated.  The  two  inverse 
transformation  matrixes  between  the  two  coordinate  systems  are  calculated 
and  stored.  The  unit  normal  n  is  used  as  the  z-axis.  The  orientation  of 
the  x  and  y  axes  is  immaterial. 

There  are  two  standard  constraints.  The  first  is  the  cylinder  of  radius  Rj 
around  the  unit  normal  n  while  the  second  is  one  of  radius  R2  also  around 
n.  Rj  =  0  indicates  that  there  is  no  inner  constraint. 

2.1.2  Cone 


A  cone  bounded  by  two  planes  perpendicular  to  the  axis,  see  Figure  2*2, 
is  used  as  the  basic  cone  configuration.  It  is  described  by  the  coordinates 
of  the  two  center  points  (X\,  Y\,  Zj),  (X2,  Y2»  Z2)  and  the  radii  Rj  and  R2. 
The  parameters  in  this  sequence  are  stored  as  PTLCTR  of  the  surface. 

In  subroutine  CONE  coefficients  of  the  surface  and  the  constraint  equations 
in  terms  of  the  system  coordinates  X-Y-Z  are  calculated.  Two  transforma¬ 
tion  matrixes  between  the  system  and  the  local  coordinates  are  calculated 
and  stored.  Note  that  the  local  coordinate  origin  is  at  the  first  center  point 
(X 1 .  Yj,  Zj)  and  the  z-axis  points  toward  the  second  center  point  (X2,  Y2. 
Z2). 

There  are  two  standard  constraints.  The  first  is  the  perpendicular  plane 
through  (X).  Y],  Z j )  while  the  second  is  the  one  through  (X2,  Y2 ,  Z2). 


*  X-Y-Z  are  the  system  coordinates  while  x-y-z  are  the  local  coordinates 
except  as  otherwise  noted. 
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2.1.3  Sphere 

A  sphere  is  best  described  by  the  center  (Xj,  Yj,  Zj)  and  a  point  on  the 
surface  (X2.  Y2,  Z2),  see  Figure  2-3.  Similar  to  a  plane,  a  sphere  is 
bounded  by  two  cylinders  of  radii  Rj  and  R2.  The  parameters  in  this 
sequence  are  stored  as  PTLCTR  of  the  surface. 

In  subroutine  SPHERE  coefficients  of  the  surface  and  the  constraint  equa¬ 
tions  in  terms  of  the  system  coordinates  are  calculated.  As  before  two 
transformation  matrixes  are  obtained.  The  origin  of  the  local  coordinates 
is  at  (Xj,  Yj,  Zj)  and  the  z-axis  points  toward  (X2,  Y2,  Z2). 

There  are  two  standard  constraints:  (1)  the  outer  cylinder  of  radius  Rj 
and  (2)  the  inner  one  of  radius  R2.  Again  R2  =  0  indicates  that  there  is  no 
inner  constraint. 

2.1.4  Paraboloid 


A  paraboloid  is  described  by  two  points  on  the  axis  of  revolution,  see 
Figure  2-4.  The  first  point  is  the  vertex  (Xj,  Yj,  Z\)  while  the  second 
is  the  focal  point  (X2.  Y2.  Z2).  Similar  to  the  plane  and  the  sphere  con¬ 
figurations,  two  standard  cylindrical  constraints  are  obtained  from  two 
radii  Rj  and  R2.  The  parameters  are  stored  in  this  sequence  as  PTLCTR 
of  the  surface. 

In  subroutine  PARABL  coefficients  are  obtained  for  the  surface  and  the 
constraint  equations.  Two  transformation  matrixes  between  the  system 
and  the  local  coordinates  are  calculated.  The  origin  of  the  local  coordin¬ 
ates  falls  on  the  vertex  (Xj,  Yj,  Z\)  and  the  z-axis  points  toward  the  focal 
point  (X2.  Y2,  Z2). 

There  are  two  standard  constraints  given  by  the  cylinders  of  radii  Rj  and 
R2.  Again  R2  =  0  indicates  that  there  is  no  inner  constraint. 

2.1.5  Plane  Baffle 

A  thin  plate  with  sharp  edged  circular  aperture  can  be  represented  by  four 
surfaces  as  shown  in  Figure  2-5.  The  first  surface  is  a  hyperboloid  of 
half  width  6.  Besides  the  radius  of  the  aperture  R,  the  surface  equation 
can  be  chosen  to  satisfy  the  edge  radius.  This  is  important  because  it  is 
the  sharpness  of  the  edge  that  determines  the  amount  of  reflection  off  a 
knife  edge  aperture. 
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FIGURE  2-3.  SPHERE  CONFIGURATION 
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The  second  surface  is  a  cone  representing  the  blade  of  the  knife  edge.  The 
remaining  two  plane  surfaces  are  two  side  faces  of  the  baffle.  The  values 
of  the  edge  radius  Re,  the  hyperboloid  half-width  6,  the  blade  angle  and  the 
baffle  thickness  are  read  in  as  a  set  of  constants. 

In  addition  to  the  above  constants,  the  parameters  necessary  to  describe 
a  baffle  are  similar  to  those  of  a  plane.  The  first  three  values  of  PTLCTR 
refer  to  the  center  of  the  aperture  (Xc,  Yc,  Zc),  and  the  second  refer  to 
the  unit  normal  of  the  aperture  (a,  (3,  y)  (see  Figure  2-5).  Note  that  the 
unit  normal  is  on  the  opposite  side  of  the  blade  surface  2.  The  radius  of 
the  aperture  R  and  the  angle  of  inclination,  p  in  Figure  2-6,  which  is  zero 
for  a  plane  baffle,  are  the  last  two  parameters  stored  as  PTLCTR  of  the 
baffle.  Since  there  are  four  surfaces  referred  to  by  a  set  of  PTLCTR,  it 
is  necessary  to  read  in  the  parameter  for  the  first  surface  only. 

In  subroutine  BAFFLE,  coefficients  are  calculated  for  the  surface  and  the 
constraint  equations.  There  are  two  standard  constraints  for  eacn  of  the 
first  two  surfaces.  Each  of  the  remaining  surfaces  has  only  one.  It  is 
sometimes  necessary  to  describe  an  outer  constraint  for  each  of  them  in 
terms  of  special  constraints. 

2.1.6  Conical  Baffle 


The  parameters  of  a  conical  baffle  are  arranged  in  exactly  the  same  way 
as  that  of  a  plane  baffle.  A  non -zero  value  is  given  for  the  angle  of 
inclination,  p  (see  Figure  2-6). 

Instead  of  hyperboloid,  the  first  surface  is  now  represented  by  a  toroid. 
This  gives  a  better  representation  of  the  knife  edge  quality,  although  it 
results  in  slower  execution  due  to  its  more  complicated  equation.  The 
second  is  again  a  cone  representing  the  blade.  The  last  two  surfaces  are 
cones  instead  of  planes. 

Again  there  are  two  standard  constraints  for  each  of  the  first  two  surfaces, 
but  only  one  for  each  of  the  remaining  two.  If  necessary,  additional  outer 
constraints  are  specified  as  special  constraints  for  surfaces  3  and  4. 


2.2  CONSTRAINTS 

The  boundary  of  a  surface  is  defined  by  a  set  of  constraints.  Figure  2-7 
shows  a  surface  with  one  of  its  constraint  surfaces.  In  order  to  test 
whether  an  intercepting  point  P  (Xj,  Yj,  Zj)  is  within  the  boundary,  we 
substitute  (X j ,  Yj,  Zj)  in  the  equation  that  defines  the  constraint  surface. 
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FIGURK  2-6.  CONICAL  BAFFLE 
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FIGURE  2-7.  CONSTRAINT 


F  (X,  V,  Z).  The  sign  of  this  value  indicates  on  which  side  of  the  constraint 
surface  the  point  P  is.  It  is  necessary  to  predetermine  and  store  the  sign 
of  P  in  the  region.  A  region  is  the  side  of  a  constraint  where  the  system 
surface  exists. 


There  are  three  types  of  constraints  as  follows: 


Linear 

Quadratic 

Angular 


Plane 

Cone,  sphere,  paraboloid,  hyperboloid,  etc. 
Toroid 


The  angular  constraint  is  specified  in  subroutine  BAFFLE  and  used  in 
subroutine  GEOM.  Since  no  input  data  is  required,  we  will  bypass  its 
discussion.  The  other  two  types  of  constraints  arc  presented  below 
followed  by  a  detailed  discussion  of  special  constraints. 
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The  type  of  a  constraint  surface  and  the  sign  of  the  constraint  function 
value  in  the  region,  for  the  Ith  constraint  of  surface  J,  are  designated  by 
the  parameter  KSTRT(I,J)  as  follows: 

KSTRT  Type  of  Constraint  and  Sign 

1  Linear  constraint  with  F  <  0  in  the  region 

2  Linear  constraint  with  F  >  0  in  the  region 

3  Quadratic  constraint  with  F  <  0  in  the  region 

4  Quadratic  constraint  with  F  >  0  in  the  region 

2.2.1  Linear  Constraint  -  Plane 

Referring  to  the  Subsection  2,1.1,  we  know  that  a  plane  is  described  by  a 
point  on  the  surface  and  a  unit  normal.  In  subroutine  PLANE  the  coefficients 
of  the  plane  are  chosen  such  that  F  >  0  on  the  side  where  the  unit  normal  is 
pointing.  Thus,  in  using  an  existing  plane  as  a  constraint,  assign 

KSTRT  =  1  when  the  unit  normal  of  the  plane  is  pointing  away  from 
the  region 

2  when  it  is  pointing  into  the  region. 

2.2.2  Quadratic  Constraints 

A  general  quadratic  constraint  is  represented  in  the  system  coordinates 
by  the  equation: 

F(X,Y,Z)  «  CjX2+C  Y2+C3Z2  +  C4XY+C5XZ+C6YZ+C7X 

+  C8Y+C9Z+C.0'°  <2'2> 

which  is  obtained  from  transforming  the  equation  in  the  local  coordinates, 

ffx.y.z)  .  c , x2  +  Cgy2  +  c3*2  +  c^xy  +  c&xz  +  c^yz  +  c?x  +C((y  +  cgz  +  c  10  =  0 

Because  of  the  second  degree  terms  in  the  equations,  the  sign  of  the  func¬ 
tion  value  is  invariant  to  transformation.  In  other  words,  for  any  point  in 
space  P  (X  j .  Y  ,  Z  )  «  p  (x  ,  y  ,  z  J )  we  have 

Sign  of  F  (Xj,  Y. ,  Zj)  »  Sign  of  f  (Xj,  y1§  Zj), 
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This  makes  it  very  easy  to  identify  the  orientation  of  a  quadraMc  surface 
by  the  following  arrangement. 

a.  Cone ,  In  subroutine  CONE,  we  use  the  z-axis  as  the  axis  of  the  cone  (see 
Figure  2-2).  Thus  both  c.j  and  C2  are  positive,  and  if  we  substitute  the 
coordinates  of  any  point  on  the  z-axis  (0,  0,  z)  in  the  equation,  we  will 
obtain 

f  (0,  0,  z) <  0. 

Therefore,  we  have 

f  <  0  and  F  <  0  for  any  point  inside  the  cone 

and 

f  >  0  and  F  >  0  for  any  point  outside  the  cone. 

b.  Sphere.  In  subroutine  SPHERE,  we  start  with  the  equation 

2  2  2  2 

f  (x,  y,  z)  =  x  +  y  +  z  -  R  =0 
Thus  we  have 

f  <  0  and  F  <  0  for  any  point  inside  the  sphere, 
f  >  0  and  F  >  0  for  any  point  outside  the  sphere. 

c.  Paraboloid.  The  paraboloid  equation  is  written  as 

f  (*•  z)  ’  *  +y2  +  v  * 0 

in  subroutine  PARABL  (see  Figure  2-4).  By  keeping  a  constant  z  and 
decreasing  the  magnitude  of  x  and  y,  we  obtain 

f  <  0  and  F  <  0  for  any  point  "inside"  the  paraboloid  (see  Figure  2-8) 

and  consequently. 

f  >  0  and  F  >  0  for  any  point  "outside"  the  paraboloid. 


2-11 


( 


2-12 


T 


336G-008 


d.  Hyperboloid.  Referring  to  Figure  2-5,  we  see  that  the  hyperboloid 
surface  is  in  the  form 

t(x,  y,  7.)-.  c/  *y2.y2u1(  =  0 

whe  re 


>  0 


and 


c3<  0 

By  holding  z  constant  and  decreasing  the  magnitude  of  x  and  y,  we  obtain 

f  <  0  and  F  <  0  for  any  point  in  the  region  including  the  center  (see 
Figure  2-8), 


f  >  0  and  F  >  0  for  any  point  in  the  opposite  region. 

The  above  observation  is  summed  up,  together  with  that  of  a  plane,  in 
Figure  2-8.  The  unit  normals  are  included  to  assist  in  memorizing  the 
orientation  of  the  surfaces.  They  always  point  toward  the  positive  sides 
of  the  surfaces. 

Modification  of  the  set  of  standard  constraints  given  by  the  basic  configu¬ 
rations  (see  Subsection  2.1)  is  accomplished  by  specifying  the  arrays 
NSTRT,  KSIDSP,  KSTRT,  PTLCSP  and  CSTRT  in  Namelist  NAME3. 

NSTRT(I)  is  the  number  of  constraints  for  surface  I.  The  default  value  is 
given  by  the  basic  parameters  of  the  surface  PTLCTR  (sec  Subsection  2. 1). 

KSIDSP(J,I)  specifies  the  Jth  constraint  of  surface  I  in  the  following 
manner: 

KSIDSP  Constraint  Surface 

-1  Given  by  array  PTLSCP 

0  Standard 

Positive  Surface  number  given  by  KSIDSP 
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When  the  standard  constraint  is  to  be  left  intact,  the  parameter  KSIDSP 
assumes  the  value  zero  as  given  by  default.  When  the  constraint  is  to  be 
replaced  by  an  existing  surface,  the  number  of  the  surface  to  be  used  as 
the  new  constraint  is  specified  by  KSIDSP.  When  no  existing  surface  is 
available  for  the  constraint,  KSIDSP  is  set  at  -1. 

KSTRT(J.I)  is  the  type  of  the  Jth  constraint  of  surface  I.  When  the  type 
of  a  constraint  is  different  from  the  standard  type  (see  Subsection  2.1), 
the  parameter  KSTRT  is  specified  as  follows: 

Existing  Surface: 

KSTRT  =  1  for  plane  with  normal  pointing  away  from  the  region 

2  for  plane  with  normal  pointing  into  the  region 

3  for  quadratic  with  normal  pointing  away  from  the  region 

4  for  quadratic  with  normal  pointing  into  the  region 

New  Surface: 

KSTRT  =  1  or  2  for  plane 

3  for  quadratic  with  normal  pointing  away  from  the  region 

4  for  quadratic  with  normal  pointing  into  the  region 

The  parameters  of  a  new  surface  as  the  Jth  constraint  of  surface  I  are 
specified  by 

PTLCSP(K.J.I)  K  =  1,  2,  ...,  6  or  8 

The  parameters  have  the  same  meaning  as  those  of  PTLCTR  for  basic  con¬ 
figurations  described  in  Subsection  2.  1.  However,  there  are  only  six  param¬ 
eters  instead  of  eight  for  a  plane,  a  sphere  or  a  paraboloid  because  no  radii 
are  raquired  for  the  constraints. 

The  above  can  be  best  demonstrated  by  the  following  example.  In  Figure 
2-9  we  have  surface  17  to  be  bounded  by  the  following  four  constraints -- 
existing  cone  of  surface  number  34,  existing  paraboloid  of  surface  number 
26,  existing  plane  of  surface  number  43,  and  a  new  surface. 

Referring  to  Figure  2-8,  we  obtain  the  directions  of  the  unit  normals  for 
the  three  existing  surfaces.  The  types  of  constraints  are  determined  as  3, 

4  and  1,  respectively,  for  these  surfaces. 
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CONSTRAINT  EXISTING  SURFACE  TYPE  OF 

NUMBER  NUMBER  CONSTRAINT 

1  34  3 

2  26  4 

3  43  1 

4  -1  I 


FIGURE  2-9.  SPECIAL  CONSTRAINT  EXAMPLE 
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The  new  surface  is  determined  by  a  point  on  the  surface  (1.2,  2.5,  3.2) 
and  the  unit  normal  pointing  into  the  region  (-0.6,  0.8,  0).  The  input  data 
thus  appears  as  follows: 

NSTRT  (17)  =  4 

K  SI  DSP  (1,  17)  =  34.  26,  43,  -1 
KSTRT  (1,  17)  =  3,  4,  1,  1 

PTLCS1M1.  4,  17)=  1.2,  2.5.  3.2,  -0.6,  0.8,  0 

When  a  constraint  surface  is  required  for  which  there  is  no  available  basic 
surface,  the  coefficients  of  the  constraint  equation.  Equation  2-2,  can  be 
read  in  through  the  array 

CSTRT(K,J,I)  =  CK  K  =  1 ,  2 . 10 

for  the  Jth  constraint  of  surface  I. 

In  the  development  of  the  program,  the  last  term  of  the  equation  was  placed 
on  the  right  hand  side.  The  user  should  change  the  sign  of  Cjq  when  he 
inputs  the  coefficients  of  Equation  2-2. 


2.3  SYSTEM  MODELING 

A  tremendous  gain  in  efficiency  of  the  program  is  achieved  by  dividing  an 
optical  system  into  sections.  Whether  they  are  separated  by  real  surfaces 
such  as  apertures  and  lenses  or  by  artificial  surfaces,  reduced  numbers  of 
surfaces  in  each  of  the  sections  result  in  fewer  calculations  in  search  of  an 
intercept. 

Thus,  in  describing  a  system,  we  first  divide  it  into  a  convenient  number 
of  sections.  In  each  of  the  sections,  the  surface  numbers  are  assigned 
that  relate  them  to  the  basic  surfaces  generated. 

In  addition  to  the  surface  number,  several  pertinent  parameters  are 
assigned  to  each  surface.  These  are  described  in  full  detail  in  the  sequel. 

2.3.  1  Surface  Number  -  ISFACE 


The  surface  number  refers  to  the  basic  surface  that  has  been  specified  in 
namelist  NMGEOM.  '1  he  information  related  to  the  surface  number  includes: 

a.  The  degree  of  polynomial  of  the  surface  equation  -  IDGREE. 
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b.  The  coefficients  of  the  surface  equations  -  COEF. 

c.  Number  of  constraints  -  NSTRT. 

d.  Type  of  each  of  the  constraints  -  KSTRT. 

e.  Coefficients  of  the  constraint  equations  -  CSTHT. 

f .  The  surface  radiation  coefficients  -  COAT. 

When  a  standard  basic  surface  is  used,  the  first  five  parameters  are  auto¬ 
matically  computed  by  the  related  subroutine.  When  a  set  of  constraints  other 
than  the  standard  constraints  is  required,  the  parameters  "c"  through  "e" 
are  obtained  by  prescribing  special  constraints  according  to  the  steps  given 
in  Subsection  2.2.  However,  the  first  five  parameters  can  be  read  in  as 
part  of  the  input  data. 


When  me  standard  ba^ic  suriat  e  is  a  battle,  four  consecutive  IS  FACE 
i.umoers  are  required  tor  the  tour  sui  fares  defining  the  baffle  geometry. 
The  user  inputs  only  the  first  number  and  one  PTLCTH  (see  Subsections 
3.  i.  5  and  2.  1.  6|.  The  remaining  three  numbers  are  then  automatically 
associated  with  the  baffle  and  cannot  be  used  with  other  surfaces. 

The  first  parameter  IUGHEE  is  used  to  describe  the  type  of  equation 
representing  the  surface.  The  following  three  types  of  equations  have  been 
developed  in  this  program. 

IOGHEE 

1 

2 


4 

An  IDGREE  of  value  4  is  restricted  to  use  with  a  toroid  surface  rather  than 
the  general  bi-quadratic  equation. 

1  he  surface  radiation  coefficients  are  input  in  accordance  with  the  model 
description  given  in  Subsection  3.3. 


C7*  '  C8V  *  C0Z  *  C!0 

(  *  c3y3  *  CjZ2  t  C,XY  +  c5xz  i  c6yz  ♦  C  X 

*  C8V  *  V  “  f  .Q 


Toroid  equation 
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2.  3,  2  Surface  Type  -  ISFTP 


When  a  ray  is  intercepted  by  a  surface  in  the  section,  it  has  to  be  determined 
first  whether  the  ray  is  to  be  reflected  back  into  the  same  section  or  trans¬ 
mitted  into  an  adjacent  one.  In  some  cases,  e.  ,  lens  surfaces,  the  ray 
can  be  split  in  both  directions.  The  parameter  surface  type  carries  such 
informations  in  the  following  manner: 


ISFTP 

1  Transmitting  surface,  e.  g. ,  aperture,  toward  entrance, 
i.  e.  ,  to  the  next  lower  numbered  section 

2  Transmitting  surface  toward  exit,  i.  e. ,  .o  the  next  higher 
numbered  section 

3  Transmitting  and  reflecting  surface,  e.  g. ,  lens,  toward 
entrance 

4  Transmitting  and  reflecting  surface  toward  exit 

5  Reflecting  surface,  i.  e. ,  non-transmitting 


The  section  numbers  are  assigned  in  the  ascending  sequence  from  the 
entrance  to  the  exit.  In  other  words,  the  first  section  always  has  the 
entrance  aperture  as  one  of  its  surfaces.  Its  adjacent  section  is  Section  2 
and  so  on  until  the  last  section  where  the  exit  aperture  or  a  detector  is 
located. 


When  the  intercepting  surface  has  a  surface  type  1  or  3  and  the  ray  is  to 
be  transmitted,  the  ray  is  directed  into  the  next  section  toward  the  entrance. 
Likewise,  the  ray  is  directed  into  the  next  section  toward  the  exit  when  the 
surface  type  is  2  or  4. 

2.3.3  X  Coefficient  -  KCQEF 


An  intercept  between  a  ray  and  a  surface  is  obtained  by  solving  the  simul¬ 
taneous  equations  of  the  ray  and  the  surface.  Theoretically,  the  point  of 
interception  thus  obtained  is  on  the  surface.  Unfortunately,  the  set  of 
numbers  available  to  a  computer  is  only  a  finite  subset  of  the  entire  set  of 
real  numbers  and  precise  solution  is  exception  rather  than  rule. 

Frequently  the  calculated  point  of  interception  falls  a  short  distant  in  front 
of  or  behind  the  surface.  No  particular  problem  is  caused  by  the  first 
situation.  However  when  the  point  is  behind  the  surface,  the  ray  can  be 
forced  to  reflect  off  the  surface  twice  instead  of  once. 
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This  is  best  illustrated  by  Figure  2-10.  A  ray  at  point  P  is  pointing  in  the 
direction  p.  Looking  for  interception  of  this  ray  by  the  surface  S,  we  solve 
the  pertinent  simultaneous  equations.  Due  to  round-off  errors  in  the  cal¬ 
culation,  the  intercepting  point  fall  a  short  distance  behind  S. 

The  command  then  goes  on  to  obtain  the  direction  q  of  the  reflected  ray 
and  searches  for  the  subsequent  interception.  As  is  clearly  seen  in  the 
figure,  the  ray  is  intercepted  immediately  by  the  same  surface  at  point  R. 
After  another  reflection  at  this  point,  the  ray  in  effect  penetrates  the  sur¬ 
face  even  though  S  is  a  non-transmitting  surface. 

This  is  only  one  of  many  examples  that  a  ray  is  forced  into  a  wrong  path 
due  tu  round  off  error  in  the  calculations.  However,  we  will  refrain  from 
exhausting  all  the  other  possibilities  and  be  content  with  describing  a 
practical  remedy  which  has  been  mechanized. 

In  the  Cartessian  coordinates,  we  represent  a  surface  by  the  general 
equation: 

F  (X,  Y,  Z)  =  0 

If  we  substitute  the  coordinates  of  a  point  Pi(Xj,  Yj,  Zj)  into  the  equation 
and  obtain 

F  (Xj ,  Yj,  Z{)  ±  0 

we  know  that  Pj  is  not  on  the  surface. 

Furthermore,  we  can  show  mathematically  that  in  the  vicinity  of  the  seg¬ 
ment  of  surface  S,  F  has  a  constant  sign  on  each  side  of  the  surface.  In 
other  words,  if  F  (Xj,  Yj,  Zj)<  0,  F<  0  for  all  points  on  this  side  of  the 
surface  S  and  F  >  0  for  all  points  on  the  opposite  side. 

This  is  where  the  K  coefficient  comes  into  play.  The  coefficient  is  defined 
as  follows: 

KCOEF 

0  Function  value  of  the  surface  equation  is  negative  in 

the  region. 

1  Function  value  of  the  surface  equation  is  positive  in 

the  region 
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FIGURE  2-10.  DOUBLE  REFLECTION  DUE  TO  ROUND-OFF  ERROR 
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A  region  is  used  here  to  represent  the  spatial  domain  that  iu  inside  the 
section. 

In  the  above  example,  we  have 
F  tXr  Yr  Zj)  <  0 

where  I’j  is  inside  the  region;  therefore. 

KCOEF  *  0. 

When  the  program  finds  that  point  It  is  the  no  t  intercept  and  that  it  is  at  a 
rather  short  distance  from  point  Q,  it  tests  the  validity  of  this  interception. 
The  coordinates  of  point  Q  are  substituted  into  the  function  F.  The  sign  of 
this  value  compared  with  KFOEF  will  immediately  tell  whether  point  Q  is 
inside  the  region  or  not.  In  this  particular  example,  Q  is  outside  the 
region,  and  the  interception  at  point  R  is  ignored. 

2  .  3.  A  Importance  Sampling  Surface  -  IMPSF 

While  the  importance  sampling  technique  will  be  discuBscd  in  full  detail  in 
Section  4,  the  procedure  of  assigning  importance  sampling  surfaces,  which 
we  will  simply  call  important  surfaces,  will  be  briefly  described  here. 

In  each  section,  there  can  be  up  to  three  important  surfaces.  The  param¬ 
eter  IMPSF'  given  for  each  surface  indicates  whethei  the  rays  from  the 
surface  will  be  split  such  that  one  or  more  rays  will  be  directed  towards 
the  important  surfaces.  Such  a  decision  is  specified  in  the  following 
manner. 

IMPSF 

0  No  importance  sampling 

1  Importance  sampling  toward  important  surface  1 

2  Importance  sampling  toward  important  surface  2 

3  Importance  sampling  toward  important  surface  3 

•l  Importance  sampling  toward  important  surfaces  1  or  2 

5  Importance  sampling  toward  important  surfaces  1  or  3 

6  Importance  sampling  toward  important  surfaces  2  or  3 

7  Importance  sampling  toward  important  surfaces  1,  2,  or  3 

Although  more  detailed  discussion  will  be  presented  later,  it  should  be 
pointed  out  here  that  the  above  splittings  are  in  addition  to  those  that  cover 
specular  components  of  the  radiation  coefficients. 
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When  IMPSF  »  0  for  the  surface,  the  ray  will  not  be  split.  When 
IMPSF  ■  1,  2  or  3,  the  ray  is  split  toward  the  particular  important  surface 
named.  When  IMPS/  ■  4,  5,  6  or  7,  a  uniform  random  number  is  used  to 
decide  which  of  the  important  surfaces  the  splitted  important  ray  is  aimed 
at. 
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Section  3 

STATISTICAL  MODELING  OK  RADIATION 


Probability  density  functions  giving  the  probable  paths  of  the  radiant 
energy  are  the  heart  of  the  Monte  Carlo  ray  tracing  method.  Whether  it 
be  a  collimated  laser  beam,  sl.ghtly  divergent  solar  energy,  a  reflected 
beam  off  a  surface,  or  an  emission  off  a  black  body,  a  well-defined  dis¬ 
tribution  function  has  to  be  determined  before  we  can  proceed  to  obtain 
tin*  probable  path  of  the  rays. 

We  present  here  a  number  of  radiation  models  that  have  been  used  in  a 
variety  of  optical  systems.  Kowevcr,  it  should  be  emphasized  that 
uiditional  models  will  be  needed  for  certain  new  problems.  The  computer 
program  is  arranged  in  such  a  way  that  it  is  fairly  easy  to  include  additional 
models. 

Available  radiation  models  up  to  this  point  include  a  collimated  or  divergent 
beam  coming  into  a  system  through  a  circular  aperture,  emission,  absorption 
(by  a  surface),  specular,  near  specular,  and  diffuse  reflection,  refracted 
specular,  and  diffuse  transmission,  and  diffraction  through  a  circular 
aperture. 


3.  I  EXTERNAL  RADIATION 

We  consider  here  collimated  or  divergent  adtant  energy  entering  the 
system  uniformly  through  a  circular  aperture.  Since  a  collimated  beam 
can  be  obtained  by  specifying  a  zero  divergence  angle  to  a  divergent  beam, 
we  will  limit  our  discussion  to  the  latter. 

A  typical  divergent  beam  is  shown  in  Figure  3-1,  passing  through  a 
circular  entrance  aperture.  The  aperture  is  placed  on  the  x-y  plane  ol 
the  local  coordinates  x-y-z,  with  the  z-axis  passing  through  the  center. 
The  local  coordinates  can  be  in  any  orientation  with  respect  to  the  system 
coordinates  X-Y-Z. 

The  direction  of  the  l>enm  as  a  whole,  which  we  shall  call  the  center  beam, 
is  defined  by  tiie  elevation  and  the  azimuthal  angles.  The  beam  is  uni¬ 
formly  spread  over  the  solid  angle  around  the  center  beam  within  the 
diNoruencc  angle. 
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it  is  important  to  realize  that  the  elevation  angle,  $Q,  and  the  azimuthal 
angle,  0O,  are  defined  with  respect  to  the  system  coordinates.  This  is 
clearly  shown  in  Figure  3-1  where  the  coordinate  system  X'-Y'-Z*  is  a 
translation  of  the  system  coordinates  to  point  P,  the  point  of  interception 
between  the  entering  beam  and  the  aperture.  4>Q  is  the  angle  between  the 
center  beam  and  the  Z’-axis.  0Q  is  the  angle  between  the  projection  of 
the  center  beam  on  the  X’ -Y'  plane  and  the  X'-axis. 

To  generate  a  ray  statistically,  we  first  locate  point  P  on  the  aperture 
with  two  uniform  random  numbers.  The  deviation  angles  and  flj  are 
obtained  with  two  additional  random  numbers.  With  these  angles  and  the 
direction  of  the  center  beam,  calculation  of  the  ray  direction  is 
straightforward. 


3.2  EMISSION 

The  emissive  power  of  a  surface  is  given  by 
e  =  <  a  T4  A 

where*  is  the  omittance,  o  the  Stefan- Boltzmann  constant, 

T  the  temperature  and  A  the  area.  If  the  emittance  and  the  temperature 
of  a  surface  are  known,  the  emission  can  be  simulated  by  generating  a 
series  of  rays  coming  off  the  surface.  The  starting  points  of  the  rays 
should  be  uniformly  distributed  over  the  area  and  the  angular  distributions 
are  governed  by  certain  emission  models.  We  shall  first  concentrate  on 
the  method  of  picking  the  starting  points  uniformly  over  an  area. 

The  location  of  a  point  on  a  surface  is  generally  defined  by  two  coordinates, 
say  £  and  r).  In  most  cases,  the  coordinates  can  be  orthogonal  and  the  area 
bounded  by  infinitesimal  lengths  of  the  coordinates  are  simply  the  product 
of  the  sides,  namely  aA  ■  A?  An.  Under  this  condition  the  point  of  uniform 
probability  is  obtained  by  uniform  selection  of  the  coordinates.  Planes 
and  spheres  are  typical  surfaces  of  this  type.  Unfortunately  a  general 
surface  does  not  always  satisfy  this  condition. 

In  addition  to  the  above  difficulty,  excessive  calculations  are  required 
when  emission  is  to  be  considered  for  a  large  number  of  surfaces.  First, 

•  surface  will  have  to  be  selected  among  all  the  surfaces  in  accordance  with  < 

their  area  ratios  The  chore  of  calculating  areas  for  all  surface  elements 
becomes  extremely  complicated  for  a  general  program  These  difficulties 
lead  to  the  idea  of  selecting  the  points  by  following  a  statistical  procedure 
and  compensating  the  probability  with  equivalent  areas.  The  method  is 
described  in  detail  as  follows 

3-3 


I 


3366-008 


Consider  an  enclosure  of  arbitrary  shape  shown  in  Figure  3-2. 


t  0  7  2  - 


3  5  M 


; 


FIGURE  3-2.  EMITTING  ENCLOSURE 

P  is  a  point  on  the  enclosure  surface,  while  O  is  anywhere  in  the  space. 

If  an  infinitesimal  enclosure  surface  dA  is  encompassed  by  an  infinitesimal 
solid  angle  dG  as  shown  in  the  figure,  the  area  is  given  by 

dA-  ^ 

I  cos 

where  S  is  the  length  of  OP  and  <t>  is  the  angle  between  OP  and  the  surface 
normal  at  P. 

We  can  directly  apply  the  above  relation  to  obtain  the  emitting  point  and 
the  equivalent  area  of  an  emitted  ray.  First,  we  place  a  point  O,  which  we 
shall  designate  as  a  pole,  anywhere  in  the  system  and  generate  a  series 
of  pilot  rays  out  of  it  uniformly  over  all  directions.  For  each  of  these 
rays,  we  search  for  interceptions  by  the  system  surfaces. 
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Each  of  these  intercepts  is  an  emitting  point  with  the  equivalent  emitting 
area  proportional  to  the  factor  S^/ 1  cos  4|  .  If  n  pilot  rays  are  initiated 
and  intercepts  are  sought  in  both  directions  of  the  pilot  rays,  the  equivalent 
area  of  an  emitting  ray  is 

A  -  2ir  s2 
e  n  |  cos  4>| 


Thus,  the  emissive  power  associated  with  the  emitted  ray  is 


2jt  S 2  c  a  T4 
n  |  cos  <j>  | 


Clearly,  the  above  expressions  are  valid  regardless  of  the  shape  of  the 
system  surfaces.  Furthermore,  there  can  be  many  or  no  intercepts  for 
each  of  the  pilot  rays. 

An  interesting  consideration  arises  in  the  selection  of  a  pole  location.  It 
is  clear  from  the  above  discussion  that  the  probability  of  a  surface  inter¬ 
cepting  a  pilot  ray  is  proportional  to  the  factor  |  cos  d|  / S^.  More  often 
than  not,  we  are  particularly  interested  in  the  emission  of  certain  surfaces 
and  would  like  to  have  more  rays  simulating  the  emission  off  these  surfaces. 
This  is  easily  achieved  by  placing  the  pole  at  short  distances  from  these 
surfaces. 

In  some  systems,  the  emitting  surfaces  are  so  far  apart  that  it  is  desirable 
to  initiate  pilot  rays  from  different  locations.  Under  these  circumstances, 
as  many  poles  as  necessary  are  placed  at  optimum  locations  and  the  system 
surfaces  are  divided  into  groups,  each  corresponding  to  a  pole. 

The  emissive  power  of  a  surface  is  assumed  diffusely  distributed  over  the 
hemisphere.  The  direction  of  a  ray  is  selected  in  the  same  way  as  that 
of  a  diffuse  reflection.  Also  the  splitting  of  rays  for  importance  sampling 
can  be  treated  in  the  same  manner  as  in  diffuse  reflection.  Discussion  of 
this  matter  is  included  in  Subsection  3.  3. 
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3.  3  ABSORPTION- REFLECTION-TRANSMISSION 


The  tracing  of  a  ray  on  its  interception  by  a  surface  is  best  treated  by 
dividing  the  probability  into  rnr  >•  .,ents;  namely  that  of  absorption,  reflection, 
and  transmission.  Absorption  can  be  easily  accounted  for  by  reducing  the 
energy  of  the  ray.  A  random  number  is  then  used  to  choose  between 
reflection  and  transmission. 

The  direction  of  the  reflected  or  transmitted  ray  is  governed  by  certain 
probability  density  functions.  Calculations  will  be  greatly  simplified  if 
we  can  express  these  functions  in  analytical  forms. 

For  most  reflecting  surfaces,  the  distribution  of  energy  can  be  approximated 
by  specular  and  diffuse  components  (Reference  3).  The  specular  component 
is  a  collimated  beam  leaving  the  surface  in  the  specular  direction.  The 
energy  of  the  diffuse  component  is  uniformly  distributed  over  the  entire 
hemisphere  and  its  probability  density  is  given  by  Lambert's  cosine  law. 

e  «  e  cos  <4 
n 

where  en  is  the  probability  density  in  the  normal  direction  and  <b  the  angle 
between  the  normal  and  the  emerging  ray. 

These  two  components  cover  the  extreme  cases  of  perfect  mirrors  and 
gray  bodies.  In  most  real  surfaces,  however,  there  is  a  concentration  of 
energy  within  a  finite  solid  angle  around  the  specular  beam  with  the  energy 
dropping  exponentially  from  its  peak  at  the  specular  direction.  The 
probability  density  of  this  near- specular  component  is 

*  -Kq 
const •  e 

where  K  is  an  angular  decay  constant  and  a  is  the  angle  between  the 
emerging  ray  and  the  specular  direction. 

While  one  near- specular  component,  together  with  the  specular  and  the 
diffuse  components,  is  sufficient  to  approximate  most  probability  density 
functions,  it  is  sometimes  necessary  to  include  two  such  components  with 
different  angular  decay  constants.  This  will  be  demonstrated  by  an  example 
later  in  this  section. 
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Transmitted  energy  through  a  surface  is  usually  well  approximated  by  a 
specular  ray  and  diffuse  scattering.  The  specular  ray  is  refracted 
according  to  Snell's  law,  which  we  will  discuss  in  more  detail  in  the  next 
section.  The  diffuse  scattering  follows  Lambert's  cosine  law  on  the  back 
side  of  the  surface. 

Adding  up  all  the  components,  we  have 
7 

C.  =  1 

l 

where 

Cj  =  absorption  coefficient 

3  specular  reflection  coefficient 
Cg  =  diffuse  reflection  coefficient 

■  first  near  specular  reflection  coefficient 
C,.  =  second  near  specular  reflection  coefficient 
Cg  =  specular  transmission  coefficient 
C ^  =  diffuse  transmission  coefficient 

The  probability  density  distribution  functions  for  the  diffuse  and  the  near- 
specular  components  are  given  by 

-K .  -K 

fp  =  A  cos  d  +  Be  a+Ce  a  (3-1) 


l 


and 


f^  =  D  cos  d'  (3-2) 

where  fr  and  ft  are,  respectively,  the  reflection  and  the  transmission  distri¬ 
bution  functions,  <}>  is  the  angle  between  the  reflected  ray  and  the  normal, 

<f>'  is  the  equivalent  of  <£  for  transmission,  a  is  the  angle  between  the  re¬ 
flected  ray  and  the  specular  ray,  and  Ki  and  K2  are  two  angular  decay 
constants.  A,  B,  C  and  D  are  constants  determined  from  the  coefficients 
Cj  in  the  following  manner: 
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The  total  energy  of  a  diffuse  reflection  is  obtained  by  integrating  the  density 
function  over  the  hemisphere.  Hence,  we  have 

2  n  ir/2 

C 3  =  /  /  (A  cos  d)  sin  <i>  d<$  d0 

o  o 

where  0  is  the  azimuthal  angle.  Integration  gives 


V  *•  A 


(3-3) 


Similarly,  we  have 

2j r  *1 2  -K 

C4  “  /  /  <B  e  > sin  a  da  d6 , 

o  o 

a  good  approximation  as  long  as  e  is  negligible  when  a  is  replaced  by 
the  minimum  angle  between  the  specular  ray  and  the  surface  tangent. 

Integration  gives 


C4 


-K  W2 

2  7T  B(  1  -  K1  e  ) 

1  +  K.2 


For  Ki  »  5,  this  can  be  approximated  by 


'4  2 

1+Kj 


(3-4) 


Similarly  for  the  second  near-specular  scattering,  we  have 


'5  2 

1+K2 


(3-5) 
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A  slight  modification  is  introduced  here  for  the  near-specular  component. 

In  the  model  presented  above,  the  radiant  density  drops  down  exponentially 
from  its  peak  at  the  specular  direction.  In  most  cases,  the  angular  decay 
constant  is  so  large  that  the  density  is  negligible  outside  a  fairly  small  a 
angle.  However,  when  the  angular  decay  constant  is  small  and  the  incident 
ray  is  close  to  the  surface  tangent,  a  significant  amount  of  energy  would 
go  into  the  surface  according  to  the  density  function.  To  prevent  this 
situation  from  occurring,  the  probability  distribution  function  is  normalized 
by  one  of  two  methods  selected  by  the  user  through  the  parameter  IEXP. 

Referring  to  Figure  3-3,  we  have  a  ray  incident  on  a  surface  at  point  O. 

OP  is  the  specular  direction  and  plane  jr  is  perpendicular  to  it  at  point  P. 
OQ  is  the  direction  under  consideration  and  a  is,  as  defined  previously, 
the  angle  between  OP  and  OQ.  The  plane  formed  by  OPQ  is  intersected 
by  the  tangent  plane  at  OR.  The  angle  <  POR  is  the  maximum  value  of 
the  angle  a  in  that  direction  before  the  ray  will  enter  the  surface. 

If  IEXP  =  O,  the  probability  distribution  function  for  a  is  normalized  by 

P(a)=  J  C  Be  "Kla  +  Ce  'K2a  sin  a  da 

J  o  a  max 


where  C 


a  max 


/ 


amax  -K ,  a  ^  _  -K2  o  .  .  , 

(Be  1  +  Ce  ^  )sm  ada 


and  Equation  (3-1)  becomes 

fr  =  A  cos  4  +  C  Be  +  Ce 

r  a  max 


(3-6) 


Since  the  exponential  factor  is  independent  of  the  direction  of  PQ,  this 
is  referred  to  as  "symmetric"  near-specular  scatter. 


If  IEXP  =  1,  a  new  variable  given  by 


€ 


JTQ 


2a 

max 


is  substituted  for  a  in  Equation  3-1  to  give 

-K  S  -K  ? 

f  *  A  cos  <f>  +  B  e  +  C  e 
r 


(3-7) 
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FIGURE  3-3.  MODIFIED  NEAR-SPECULAR  REFLECTION  MODEL 
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The  factor  n  / 2  is  included  here  so  that  4  =  a  for  normal  incidence. 

Since  the  exponential  factor  is  dependent  on  the  direction  PQ  through  the 
variable  4,  this  is  referred  to  as  "skewed'1  near-specular  scatter. 

The  input  parameters  for  a  surface  are  summarized  in  Table  3-1. 

The  following  example  of  a  typical  high  quality  mirror  is  presented  to 
demonstrate  the  usage  of  the  various  components. 

Example:  Instead  of  attempting  to  match  a  measured  angular  radiant 
energy  distribution  by  trial  and  error,  we  shall  assume  a  set  of  parameters 
and  plot  the  resultant  distribution  of  a  typical  high  quality  mirror.  Using  the 
symmetric  near-specular  model,  we  assume  the  following: 


-6 

7T  X  10 

-6 

10 

Kj  =  500 

-7 

5  x  10 

K2  ■  250 

TABLE  3-1.  SURFACE  RADIANT  DISTRIBUTION  PARAMETERS 


Input  Data 
COAT  (I,  J)  * 

COAT  (14,  J) 
COAT  (15,  J) 
IEXP  *  0 
IEXP  ■  1 


Parameters  and  Choice 
_ of  Models _ 

£  *  J  is  the  surface  number. 

1  I  »  1.  2,  ....  7 

K1 

K2 

Symmetric  near- specular  model  given  by  Equation  3-6. 
Skewed  near-specular  model  given  by  Equation  3-7. 


Substituting  these  parameters  in  Equations  3-3,  3-4,  3-5,  and  3-6,  we 
obtain  the  individual  and  the  combined  angular  radiant  density  distributions 
plotted  in  Figure  3-4. 
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FIGURE  3-4.  RADIANT  DENSITY  DISTRIBUTION  OF  DIFFUSE  AND 

NEAR-SPECULAR  REFLECTION 
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In  the  simulation  of  a  mirror  with  a  measured  biangular  radiant  energy 
distribution,  a  different  set  of  component  coefficients  is  assumed  and 
substituted  in  the  equations  until  the  combined  density  distribution  will 
represent  the  measured  distribution  curve. 


3.4  REFRACTION 

Given  the  indices  of  refraction  for  both  sides  of  a  surface,  the  direction 
of  the  refracted  ray  is  determined  by  Snell's  law,  which  states,  "the 
refracted  ray  lies  in  the  plane  of  incidence,  and  the  sine  of  the  angle  of 
refraction  bears  a  constant  rat  o  ,c  the  sine  of  the  angle  of  incidence" 
(Reference  4).  In  the  refraction  at  a  boundary  between  two  substances 
having  indices  of  refraction  n  and  n',  it  may  be  written  in  the  form 

n  Sin  d  *  n'  Sin  d' 

where  d  is  the  angle  between  the  incident  ray  and  the  inward  normal  and 
d'  the  angle  between  the  refracted  ray  and  the  outward  normal. 

Tracing  of  a  refracted  ray  accordingly  is  straightforward,  except  in  the 
case  of  total  internal  reflection,  where 

£sin  d*  2  1 

and  the  refracted  ray  is  reflected  instead.  Lacking  a  reliable  prediction 
of  the  probable  path  of  the  ray,  we  assume  that  the  totally  reflected  ray 
is  distributed  according  to  the  reflection  components. 


3.5  DIFFRACTION 

Diffraction  of  a  bundle  of  radiant  energy  through  a  plane  aperture  of 
arbitrary  shape  has  been  successfully  modeled  (Reference  2),  Besides 
those  given  in  the  reference,  additional  tests  of  the  model  have  oeen  made, 
which  wc  will  present  ns  examples  in  Section  9. 
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The  procedure  to  be  followed  in  determining  the  direction  of  a  diffracted 
ray  is  summarized  below- 

a.  Follow  the  ray  until  it  strikes  the  aperture  plane.  Shown  in 
Figure  3-5  is  such  an  aperture  and  P  is  the  point  of  interception. 

b.  oj  is  the  shortest  distance  from  P  to  the  aperture  edge. 

c.  Using  P  as  the  center  and  0|  as  the  mirror  radius,  fit  the  biggest 
ellipse  possible  inside  the  aperture.  02  is  the  major  radius.  £  -  n 
are  coordinate  axes  in  the  directions  as  shown. 


d.  Two  angle  constants  are  calculated  from 


tan 


-1 


I 

2kO| 


and 


■  tan 


-I 


1 


where  k  <s  the  wave  number. 
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FIGURE  3-5.  TYPICAL  APERTURE  PLANE 


e.  Illustrated  in  Figure  3-6,  ftj  is  the  deviation  angle  of  the  ray  on 
the  plane  given  by  the  incident  ray  and  the  axis  P& 

Similarly.  0g  is  the  deviation  angle  in  the  Pn  direction. 

f.  The  direction  of  the  diffracted  ray  is  determined  by  the  following 
probability  density  functions. 


1  f  1  /  ^2  \  ^  1 
w— —  'xp-i(— ) 

v2i02  L  V  2  /  J 


(3-8) 


(3-9) 
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FIGURE  3-6.  DEVIATION  ANGLE  OF  RAY 


Bending  of  rays  according  to  the  above  procedure  is  available  at  circular 
apertures.  In  order  to  include  diffraction  for  certain  of  those  apertures, 
the  users  simply  specify 

KDIF  (I)  *1  for  I  *  surface  numbers. 

From  the  above  equations,  we  sec  that  the  angle  constant  6*  is  significant 
only  when  the  distance  from  point  P  to  the  aperture  edge,  o,  is  in  the  order 
of  the  wave  length  and  only  in  this  situation  that  diffraction  needs  to  be 

considered.  When  a  ray  strikes  an  aperture  plane  where  diffraction  is  | 

called  for,  the  shortest  distance  to  the  edge  a\  is  compared  with  the  wave 
length.  If  the  ratio  of  oi  to  the  wave  length  is  bigger  than  a  constant  I 

DEDGE,  specified  by  the  users,  the  diffraction  of  the  ray  is  neglected.  i 

Since  the  probability  density  is  distributed  over  a  finite  solid  angle,  a  * 

diffracted  ray  can  be  split  for  importance  sampling.  The  procedure  of  pre¬ 
scribing  such  splitting  is  in  the  general  discussion  on  importance  sampling  i 

in  the  next  section. 

J 
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Section  4 

IMPORTANCE  SAMPLING 


The  Monte  Carlo  method,  which  is  the  basis  of  the  GUERAP  simulation 
program,  is  a  statistically  accurate,  but  time-consuming  approach  to 
system  analysis.  Fortunately,  the  execution  time  can  be  reduced  several 
orders  of  magnitude  using  the  importance  sampling  technique  (Reference  1). 

When  a  ray  strikes  a  surface  where  the  energy  is  to  be  distributed  among 
several  components,  as  described  in  Subsection  3.3,  one  of  the  compo¬ 
nents  is  statistically  chosen.  If  the  energy  of  this  component  is  scattered 
over  a  finite  solid  angle,  the  direction  of  the  scattered  ray  is  randomly 
selected  according  to  the  appropriate  probability  density  function. 

In  certain  surface  radiation  models,  some  components  have  such  pre¬ 
dominant  probability  that  an  extremely  large  number  of  rays  will  have  to 
be  generated  before  any  rays  will  select  any  of  the  other  components. 

As  the  first  step  in  importance  sampling,  a  ray  is  split  into  several 
rays,  each  corresponding  to  one  component.  To  keep  the  number  of  rays 
under  control,  this  splitting  is  done  after  the  energy  is  reduced  by  absorp¬ 
tion  and  the  decision  has  been  made  between  reflection  and  transmission. 
The  splitting  is  automatically  performed  by  the  program. 

Further  importance  sampling,  which  requires  specification  by  the  user, 
is  explained  in  more  detail  in  the  following. 


4.  1  IMPORTANCE  SAMPLING  AND  SCATTERING 

Within  a  particular  component  where  the  energy  is  scattered  over  a 
finite  solid  angle,  the  direction  of  the  scattered  ray  is  randomly  selected 
according  to  the  probability  density  function.  In  most  analyses,  however, 
we  are  most  interested  in  the  amount  of  radiant  energy  reaching  certain 
surfaces,  if  the  rays  arc  allowed  to  select  their  paths  randomly,  a 
great  portion  of  the  execution  time  will  be  wasted  in  tracing  rays 
through  areas  of  little  interest  to  the  user. 

The  situation  can  be  dramatically  improved  by  the  use  of  importance  sampling. 
Similar  to  the  component  splitting  given  above,  a  non-specular  component  is 
further  split  into  multiple  rays,  each  containing  a  portion  of  the  scattered 
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energy.  The  directions  of  all  but  one  of  the  rays,  which  we  shall 
designate  as  important  rays,  are  chosen  by  the  user  to  aim  toward  the 
areas  of  interest,  and  the  fractions  of  the  energy  to  be  assigned  to  each 
ray  are  calculated  from  their  appropriate  probability  density  functions. 

The  remaining  rays  follow  the  direction  statistically  selected  to  complement 
the  entire  probability  density  function.  Its  energy  is  equal  to  the  total 
energy  loss  that  assigned  to  the  important  rays  but  it  is  restricted  to 
directions  outside  the  solid  angle  defined  for  the  important  rays. 

The  importance  sampling  technique  as  applied  to  scattering  of  energy  is 
summarized  as  follows. 

Referring  to  Figure  4-1,  we  have  energy  radiating  from  a  surface 
element  at  0.  If  Iq  is  the  total  energy  and  E  (ft)  is  the  angular  distribution 
function  of  the  energy,  we  have 

V  //  E<n> dn 

o 

It  will  become  clear  later  that  it  is  not  necessary  to  obtain  the  function 
E  (ft)  in  the  process  of  importance  sampling. 

The  total  energy  which  is  simulated  by  a  single  ray  in  the  basic  Monte 
Carlo  method  is  to  split  into  two  rays  at  0.  One  of  the  rays,  the  important 
ray,  is  aimed  at  the  surface  Aj,  the  important  surface,  specified  by  the 
user.  Although  Ai  appears  as  part  of  a  plane,  it  could  be  any  type  of 
surface  as  far  as  the  importance  sampling  technique  is  concerned. 

If  Ai  is  defined  by  the  orthogonal  coordinates  f  -  rj  and  a  point  P  in  Aj 
is  randomly  chosen  with  respect  to  the  coordinates,  OP  can  be  used  as 
the  important  ray  with  energy  given  by 

I,  *  //  E(f,  rj)  d  f  d  r) 

*  A 


or 


ij  =  e2  (? ,  n)  A1 

where  El  (?  ,  n)  is  the  average  energy  per  unit  area  of  Aj.  If  a  large 
number  of  rays  are  traced,  the  average  energy  can  be  replaced  by  the 
energy  at  point  P  and  we  have  statistically 

Ij  ■  E  (£  ,  q)  Aj. 
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The  complementary  ray  which  is  required  to  miss  surface  Aj  will  have 
the  energy 


4.  2  IMPORTANCE  SAMPLING  AND  DIFFRACTION 

Importance  sampling  of  rays  toward  specified  areas  can  be  employed 
as  soon  as  a  probability  density  function  exists  over  a  finite  solid  angle. 

This  can  be  reflection  or  transmission  off  a  physical  surface  or  diffraction 
through  an  aperture. 

As  described  in  Subsection  3.  5,  the  direction  of  a  ray  in  passing  through 
an  aperture  is  governed  by  the  probability  density  functions  given  by 
Equations  3-8  and  3-9.  As  soon  as  a  test  shows  that  the  intersection 
point  is  within  DEDGE  wavelengths  of  the  aperture  edge,  the  ray  can  be 
split  as  in  surface  scattering. 

The  same  procedure  is  followed  in  selecting  the  directions  of  the  important 
rays.  The  fractions  of  the  energy  assigned  to  these  rays  are  determined 
by  the  probability  functions.  A  ray  is  statistically  selected  to  account 
for  the  remaining  probability.  Again,  its  energy  is  equal  to  the  total 
energy  less  that  assigned  to  the  important  rays  and  it  is  restricted  to  the 
directions  outside  the  solid  angle  defined  for  the  important  rays. 

A  special  importance  sampling  technique  has  been  implemented  at  the 
entrance  aperture.  As  we  see  from  above,  appreciable  bending  of  rays 
due  to  diffraction  occurs  only  when  the  ray  is  within  a  short  distance  of 
the  aperture  edge.  In  most  optical  systems,  the  entrance  apertures  are 
so  large  that  most  of  the  rays  generated  according  to  the  model  given  in 
Subsection  3. 1  will  be  in  the  center  region. 

In  order  to  concentrate  on  the  diffraction  effect  of  the  entrance  aperture, 
it  is  divided  into  a  narrow  outer  ring  and  the  center  circle.  The  width 
of  this  ring  is  given  by  the  constant  DRENT.  When  a  positive  value  is 
specified  for  this  constant,  the  number  of  rays  generated  at  the  entrance 
aperture  are  evenly  divided  between  the  two  areas. 

The  rays  in  the  ring  are  subsequently  split  into  important  and  comple¬ 
mentary  rays  if  important  surfaces  are  specified  for  the  entrance  aperture. 
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4.  3  IMPORTANT  SURFACE  DEFINITION 

For  each  of  the  important  surfaces  corresponding  to  a  surface  where  a 
ray  is  to  be  split,  a  point  is  randomly  selected  on  the  important  surface. 

The  important  ray  is  obtained  by  proceeding  toward  this  point  and  the 
corresponding  fraction  of  energy  is  calculated  from  the  probability 
density  function. 

The  number  of  important  surfaces  in  each  section,  with  3  as  the  upper 
limit,  is  given  by  the  array  NIMPS(I),  where  I  is  the  section  number. 

Each  important  surface  can  be  further  partitioned  as  described  in  Sub¬ 
section  4.  4. 

An  important  surface  must  be  defined  as  a  plane  bounded  by  an  outer  and 
probably  an  inner  constraint.  Since  it  is  in  exactly  the  same  form  as  a 
basic  plane  described  in  Subsection  2. 1,  the  important  surfaces  are  simply 
read  in  as  basic  planes.  These  planes  that  are  generated  for  the  system 
configuration  or  additional  planes  generated  solely  as  important  surfaces. 

After  sufficient  planes  have  been  generated  and  their  equation  coefficients 
calculated,  it  is  necessary  only  to  specify  the  surface  numbers  of  those 
surfaces  to  be  used  as  important  surfaces.  These  surface  numbers  are 
given  as  array  IMPORT  with  the  first  subscript  referring  to  the  sequential 
number  of  important  surfaces  in  the  section,  and  the  second  subscript  the 
section  number.  For  example,  IMPORT  (J,  I)  ■  K  means  that  the  basic 
plane  identified  as  surface  K  is  the  Jth  important  surface  of  section  I. 

In  addition,  the  number  of  important  surfaces  in  each  section  must  be  speci¬ 
fied  in  the  array  NIMPS  with  the  subscript  referring  to  the  section  number. 

Assignment  of  the  important  surfaces  to  the  system  surfaces  is  specified 
by  the  array  IMPSF.  Users  are  referred  to  Paragraph  2.  3. 4  for  its 
usage. 


4.4  IMPORTANT  SURFACE  PARTITIONING 

In  the  brief  description  of  the  importance  sampling  technique  given  in 
Subsection  4. 1,  we  substitute  the  energy  density  at  P  for  the  average 
value,  i.  e. , 

//e  <? ,  n)  d  I  a  n  *  e  ((  .  y  A 

Ai 

This  is  statistically  correct  and  will  provide  accurate  results  if  a  large 
number  of  rays  are  traced  through  the  system. 
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In  order  to  keep  in  check  the  number  of  rays  required  it  is  advisable  to 
reduce  the  error  attributable  to  the  above  approximation.  If  we  divide 
the  surface  Ai  into  n  area  elements  and  generate  an  important  ray  for 
each  of  them  the  approximation  will  be  in  the  form 

n 

ffE(Z,n)dkdn=  E  (4  .  i) )  A 

Aj  i  =  1  * 

and  clearly  the  individual  error  will  be  drastically  reduced. 

Applying  the  above  improvement  to  the  program,  we  divide  each  important 
surface  into  equal  areas  of  rings  and  segments,  as  illustrated  in  Figure  4-2. 
This  is  specified  by  the  user  through  the  arrays  NRING  and  NSEG.  Similar 
to  IMPORT,  there  are  two  subscripts  for  each  of  the  arrays.  The  first 
subscript  refers  to  the  sequential  number  of  important  surface  in  a  section 
while  the  second  subscript  refers  to  the  section  number.  The  choice  between 
rings  and/or  segments  is  optional:  **ings  should  be  emphasized  if  the 
probability  density  distribution  is  approximately  axi-symmetric  over  the 
important  surface  and  segments  should  be  emphasized  if  it  is  not. 

How  much  an  important  surface  should  be  partitioned  depends  very  much 
on  the  probability  density  distribution.  A  rule  of  thumb  is  that  the  density 
should  vary  no  more  than  a  factor  of  5  over  each  area  element  and  the 
energy  given  to  any  important  ray  should  be  at  least  a  factor  of  3  smaller 
than  the  total  energy. 

In  cases  where  the  density  distribution  is  fairly  uniform,  e.  g.  diffuse 
emission,  the  above  criteria  is  ascertained  if  we  keep  the  area  ratio, 
i.  e. ,  the  ratio  of  the  partioned  area  to  the  square  of  the  average  length 
of  an  important  ray  to  the  area  element,  within  the  factor  0.  2.  In  other 
cases,  e.  g. ,  diffraction,  the  above  area  ratio  should  be  kept  very  small. 

The  user  should  inspect  the  energy  levels  to  ensure  that  the  above  criteria 
are  met. 


4.  5  IMPORTANT  SURFACE  DEFINITION  GUIDE 

The  key  to  efficient  use  of  GUERAP  is  a  judicious  use  of  importance 
sampling.  Much  of  the  judgment  must  be  obtained  by  experience; 
the  purpose  of  this  section  is  to  provide  a  basis  for  avoiding  the  more 
obvious  difficulties  that  the  user  may  encounter. 
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4.5.1  Importance  Surface  Selection 

Locating  an  important  surface  for  the  direct  scattered  energy  from  each 
real  surface  is  relatively  simple.  In  addition,  the  user  must  investigate 
the  possibility  that  indirect  scattered  energy  (e.  g.  multiple  reflection 
energy)  will  be  as  significant  as  the  direct  scattered  energy.  In  this 
case  the  user  should  define  a  set  of  additional  importance  surfaces  to 
"guide"  tne  indirect  scattered  energy  to  the  principal  importance  surface. 

Figure  4-3  illustrates  this  technique.  The  single  section  in  Figure  4-3 
consists  of  two  planes,  surfaces  2  and  4,  defining  a  circular  entrance 
aperture  (surface  1)  and  a  circular  exit  aperture  (surface  3)  respectively. 

The  planes  are  bounded  by  a  concentric  cylinder,  surface  5. 

Energy  striking  surface  4  cannot  be  scattered  directly  into  the  exit  aperture, 
the  principal  important  surface.  However,  a  significant  fraction  of  the 
energy  leaving  surface  4  will  strike  surface  2  and  then  be  rescattered  from 
surface  2  through  the  exit  aperture.  Surface  2  should  therefore  be  defined 
as  an  importance  surface  for  surface  4. 

4.  5.  2  Important  Surfaces  for  Diffracted  Energy 

In  many  cases,  as  energy  is  diffracted  at  an  aperture,  what  is  significant 
is  not  the  energy  that  subsequently  strikes  a  real  surface  or  exit  aperture, 
but  is  the  energy  that  is  diffracted  into  the  field  of  view  at  the  aperture. 

In  this  case  the  user  should  define  an  imaginary  important  surface  which 
subtends  a  solid  angle  representing  the  system  field  of  view  at  the  diffraction 
zone. 

Figure  4-4  illustrates  this  situation.  Only  the  energy  which  is  diffracted 
into  the  field  of  view  will  be  focused  by  the  lens  into  the  exit  aperture. 

The  exit  aperture  cannot  be  used  as  the  importance  surface;  the  lens  will 
focus  the  resulting  important  rays  in  front  of  the  exit  aperture.  The  lens 
surfaces  cannot  be  used;  very  few  of  the  resulting  rays  will  be  in  the  field 
of  view. 

The  user  must  define  the  imaginary  surface  which  essentially  defines  the 
field  of  view  at  the  diffraction  ring.  1  his  is  illustrated  in  Figure  4-5.  Two 
cones  defining  the  field  of  view  are  extended  from  opposite  sides  of  the 
diffraction  zone,  parallel  to  the  system  optical  axis.  The  cones  are  extended 
past  their  point  of  intersection  A,  and  the  important  surface  is  located  at 
B,  with  an  outer  radial  constraint  equal  to  the  overlap  radius  BC.  This 
surface  subtends  a  solid  angle  equal  to  a  fraction  of  the  field  of  view,  which 
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is  defined  by  the  outside  limits  of  the  cones,  radius  BD.  As  AB  is  made 
larger  and  larger,  the  fraction  approaches  unity.  If  AB  is  made  too  large 
however,  the  simulation  will  suffer  due  to  computer  round-off  errors.  A 
practical  choice  is  to  make  AB  just  large  enough  that  the  fraction  is  0.90 
to  0. 95. 
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Section  5 

PROCR  A  ”  DESCRIPTION 


To  assist  the  user  in  his  understanding  of  the  GUERAP  computer  program 
and  in  his  effort  to  modify  it  when  necessary,  we  present  a  concise  descrip¬ 
tion  in  this  section. 

The  program  can  be  divided  into  seven  major  portions,  each  consisting  of 
several  closely  related  subroutines.  After  a  brief  description  of  its  main 
task,  we  list  the  subroutines  and  their  corresponding  functions. 

At  the  end  of  this  section,  skeleton  flow  charts  are  given  for  the  main  pro¬ 
gram  and  its  major  subroutines.  Together  with  the  embedded  comments 
throughout  the  program  deck,  these  should  greatly  facilitate  any  desired 
user  modifications. 


5.  1  MAIN  PROGRAM  GUERAP 

a.  Sets  the  values  of  constants  and  non-zero  default  values  of  variables. 

b.  Reads  in  surface  parameters. 

c.  Calls  subroutines  PLANE.  CONE,  SPHERE,  PARABL  and  BAFFLE 
to  generate  the  basic  surfaces  and  their  standard  constraints. 

d.  Reads  in  the  parameters  of  and  calculates  the  special  constraints. 

e.  Reads  in  and  prepares  surface  radiation  coefficients. 

f.  Reads  in  parameters  describing  the  system. 

g.  Reads  in  parameters  describing  the  set  of  important  surfaces. 

h.  Reads  in  the  parameters  of  and  performs  ray  tracing  of  external 
radiation. 

i.  Prints  out  results  thereof. 

j.  Reads  in  the  parameters  of  and  performs  ray  tracing  of  internal 
emission. 

k.  Prints  out  results  thereof. 
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5.2  BASIC  SURFACE  SUBROUTINES 

These  subroutines  generate  basic  surfaces  and  handle  transformation  of 

their  equations  between  system  and  local  coordinates. 

a.  Subroutine  PLANE 

1.  Calculates  coefficients  for  the  equations  of  planes  and  their  stand¬ 
ard  constraints. 

2.  Calculates  coefficients  for  the  equations  of  plane  special  constraints. 

3.  Calls  subroutine  EQUATN  to  transform  the  equations  into  the  sys¬ 
tem  coordinates. 

b.  Subroutine  CONE 

1.  Calculates  coefficients  for  the  equations  of  cones  and  their  stand¬ 
ard  constraints.  Note  that  circular  cylinders  are  treated  as  cones. 

2.  Calculates  coefficients  for  the  equations  of  conical  special 
constraints. 

3.  Calls  subroutine  EQUATN  to  transform  the  equations  into  system 
coordinates. 

c.  Subroutine  SPHERE 

1.  Calculates  coefficients  for  the  equations  of  spheres  and  their 
standard  constraints. 

2.  Calculates  coefficients  for  the  equations  of  spherical  special 
constraints. 

3.  Calls  subroutine  EQUATN  to  transform  the  equations  into  the 
system  coordinates. 

d.  Subroutine  PARABL 

1.  Calculates  coefficients  for  the  equations  of  paraboloids  and  their 
standard  constraints. 

2.  Calculates  coefficients  for  the  equations  of  parabolic  special 
constraints. 
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3.  Calls  subroutine  EQUATN  to  transform  the  equations  into  the 
system  coordinates. 

».  Subroutine  BAFFLE 

1.  Calculates  coefficients  for  the  equations  of  surfaces  representing 
baffles  and  their  standard  constraints. 

2.  Calls  subroutine  EQUATN  to  transform  the  equations  except  those 
of  toroids  into  the  system  coordinates. 

f.  Subroutine  EQUATN 

1.  Transforms  the  above  equations  into  system  coordinates. 

2.  Generates  and  stores  the  transformation  matrixes. 

g.  Subroutine  TRNSFM 

Transforms  the  location  and  orientation  of  a  ray  between  local  and 
system  coordinates. 

5.  3  SUBROUTINES  ON  EXTERNAL  RADIATION 

a.  Subroutine  EXTRADN 

1.  Calls  subroutine  EXT  to  generate  rays  at  entrance  aperture. 

2.  Calls  subroutine  TRACING  and  traces  the  rays. 

3.  Prints  out  results  on  external  radiation. 

b.  Subroutine  EXT 

Generates  rays  at  entrance  aperture. 

5.  4  SUBROUTINES  ON  INTERNAL  EMISSION 
a.  Subroutine  EM1SSN 

1.  Generates  pilot  rays  at  poles. 
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2.  Obtains  emitting  points  through  interceptions  between  pilot  rays 
and  related  surfaces. 

3.  Calls  subroutine  BBRFN  to  calculate  black  body  radiation  function. 

4.  Generates  emitted  rays. 

5.  Calls  subroutine  TRACING  to  trace  the  rays. 

6.  Prints  out  results  on  internal  thermal  emission, 

b.  Subroutine  BBRFN 

Calculates  black  body  radiation  function  in  terms  of  wave  bands. 

5.  5  SUBROUTINES  ON  RAY  TRACING 

This  group  handles  Monte  Carlo  ray  tracing  through  system  in  detail. 

a.  Subroutine  TRACING 

1.  Traces  rays  through  the  system. 

2.  Recalls  important  rays  from  storage  and  traces  them  through  the 
system. 

b.  Subroutine  GEOM 

1.  Calculates  points  of  interception  between  rays  and  surfaces. 

2.  Tests  the  points  in  terms  of  constraints. 

3.  Selects  the  closest  interception. 

c.  Subroutine  NORM 

Calculates  direction  cosines  of  surface  unit  normals  at  points  of 
interception. 

d.  Subroutine  TRANS 

Transforms  orientation  of  rays  from  local  coordinates  at  points  of 
interception  to  system  coordinates. 
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e.  Subroutine  REFLIN 

1.  According  to  surface  radiation  model,  obtains  directions  of  rays 
after  its  interactions  with  surfaces. 

2.  Calls  subroutine  REFRN  for  refracted  rays. 

3*  Calls  subroutine  EXPECT  to  generate  and  store  important  rays. 

f.  Subroutine  REFRN 

1.  Calculates  directions  of  refracted  rays  according  to  Snell1  s  Law. 

2.  Tests  total  internal  reflection. 

g.  Subroutine  ANGSSD 

Calculates  directions  of  rays  for  specular  reflection,  diffuse  reflec¬ 
tion,  near  specular  reflection  and  diffuse  transmission. 

h.  Subroutine  SCATTR 

Calculates  off-specular  angle  for  near-specular  reflected 
rays. 

i.  Subroutine  DIFRN 

1.  According  to  diffraction  model,  obtains  directions  of  rays  after 
they  pass  through  circular  apertures. 

2.  Calls  subroutine  DIFEXP  to  generate  and  store  important  rays. 

j.  Subroutine  DISTB 

Records  and  prints  ray  count  and  energy  maps  on  exit  apertures. 

5.6  SUBROUTINES  ON  IMPORTANCE  SAMPLING 
a.  Subroutine  EXPECT 

1.  Calls  subroutine  ISELECT  to  select  important  surfaces. 

2.  Generates  important  rays  according  to  surface  radiation  models. 

3.  Calls  subroutine  RAYSAV  to  store  important  rays. 


5-5 


3366-008 


4.  Generates  random  rays  according  to  surface  radiation  models. 

5.  Calls  subroutine  TESTSS  to  test  interceptions  of  complementary 
rays  by  important  surfaces. 

b.  Subroutine  DIF  EXP 

1.  Calls  subroutine  ISELECT  to  select  important  surfaces. 

2.  Generates  important  rays  according  to  diffraction  model. 

3.  Calls  subroutine  RAYSAVto  store  important  rays. 

4.  Generates  random  rays  according  to  diffraction  model. 

5.  Calls  subroutine  TESTSS  to  test  interceptions  of  complementary 
rays  by  important  surfaces. 

c.  Subroutine  ISELECT 
Selects  important  surfaces. 

d.  Subroutine  TESTSS 

Tests  interceptions  of  complementary  rays  by  important  surfaces. 

e.  Subroutine  RAYSAV 
Stores  important  rays. 

5.  7  MISCELLANEOUS  SUBROUTINES 

a.  Subroutine  PLOT 

Plots  radial  energy  density  distribution  on  the  exit  aperture. 

b.  Subroutine  BIQUA 

Solves  the  bi-quadratic  equation, 

C^4  ■<  C2X3  +  C3X2  +  C4X  +  C5  ■  0. 
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c.  Subroutine  CUBIC 
Solves  the  cubic  equation, 

CjX3*^2  +c3x  +  c4  -  0. 

d.  Function  ACOS 
Obtain  y  ■  cos"1  X. 

5.  8  MAIN  PROGRAM  AND  SUBROUTINE  FLOW  CHARTS 

Flow  charts  are  given  for  the  main  program  and  selected  subroutines  in 
Figures  5-1  through  5-7. 
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FIGURE  5-2.  SUBROUTINE  EXTRADN  FLOW  CHART 
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RESULT  / 
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(^return) 


FIGURE  5-3.  SUBROUTINE  EMISSN  FLOW  CHART  (1  of  2) 
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FIGURE  5-3.  SUBROUTINE  EMISSN  FLOW  CHART  (2  of  2) 
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FIGURE  5-4.  SUBROUTINE  TRACING  FLOW  CHART  (1  of  3) 


5-12 


1 


3366-008 


CALCULATE  INWARD 
SURFACE  NORMAL 
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IMPORTANT  RAYS 
RETURN  WITH 
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ENERATE  G  STORE  IMPORTANT 
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RETURN  WITH  RANDOM 
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FIGURE  5-4.  SUBROUTINE  TRACING  FLOW  CHART  (2  of  3) 
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CALL  TRNSFM 


TRANSFORM  RAY  INTO 
LOCAL  COORDINATE 
k  SYSTEM  / 


CALL  DISTB 


ADD  THE  RAY 
INTO  THE  RAY 
AND  ENERGY 
MAPS 


FIGURE  5-4.  SUBROUTINE  TRACING  FLOW  CHART  (3  of  3) 
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FIGURE  5-5.  SUBROUTINE  GEOM  FLOW  CHART  (2  of  2) 
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FIGURE  5-6.  SUBROUTINE  REFLN  FLOW  CHART  (1  of  4) 


5-17 


V 


CALL  ANGSSD 
OBTAIN  ANGLES 
DEFINING  NEAR 
SPECULAR 
REFLECTION  BEAM/ 


CALL  TRANS 


'CALCULATE  DIRECTION 
OF  NEAR  SPECULAR 
REFLECTION  BEAM 


FIGURE  5-6.  SUBROUTINE  REFLN  FLOW  CHART  (2  of  4) 
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FIGURE  5-6.  SUBROUTINE  REFLN  FLOW  CHART  (3  of  4) 
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FIGURE  5-6.  SUBROUTINE  REFLN  FLOW  CHART  (4  of  4) 
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FIGURE  5-7.  SUBROUTINE  EXPECT  FLOW  CHART 
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Section  6 
INPUT  DATA 


The  complete  input  data  variables  of  the  GUERAP  computer  program  are 
listed  in  this  section.  The  data  are  divided  into  eight  groups  arranged  in 
the  sequence  given  in  the  listing. 

A  variable  can  be  either  a  scalar  or  an  array.  A  scalar  variable  is 
designated  simply  by  the  variable  name,  while  an  array  is  indicated  by  the 
name  followed  by  the  suitable  number  of  subscripts  in  a  pair  of  parentheses. 

Each  variable  is  accompanied  by  a  brief  definition,  the  default  value  and 
the  page  number  where  the  related  description  may  be  found.  While  the 
non-zero  default  values  are  given  in  the  program,  the  user  has  to  set 
the  entire  field  length  to  zero  before  execution. 

Except  for  the  first  group,  the  data  are  input  via  namelist  statements,  a 
standard  feature  in  FORTRAN  V.  The  user  is  referred  to  any  modern 
FORTRAN  manual  for  description  of  the  usage. 

Although  the  number  of  elements  in  each  array  is  limited  by  the  assigned 
dimension,  the  user  can  readily  extend  such  limits.  For  the  current  limits 
and  the  alteration  procedure,  the  user  is  referred  to  the  section  "INSTRUC¬ 
TION  FOR  CHANGING  DIMENSIONS"  in  the  listing  of  main  program  GUERAP. 

6. 1  HEADING  AND  RANDOM  NUMBER  PILOT  SEED 

Any  number  of  cards  with  alphanumeric  characters  and  blanks  can  be 
used  as  the  problem  header.  It  is  terminated  by  a  card  with  ENDTITLE 
starting  in  the  first  column.  This  header  will  appear  in  the  beginning  of 
the  output. 

Default  Reference 

Variable  Definition  Value  Page _ 

ISEED  Pilot  number  for  the  random  Must 

(in  <Z>20  format)  number  routine.  Any  odd  be 

positive  integer  less  than  2^.  input. 
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6.2  NAMELIST  NAME1 

Miscellaneous  data  for  the  system  and  external  radiation  is  presented 


below. 

Default 

Reference 

Variable 

Definition  Value 

Page 

SUNAGD 

Elevation  angle  of  incident 
beam.  * 

0 

3-2 

SUNAZD 

Azimuthal  angle  of  incident 
beam. 

0 

3-2 

SUNDVD 

Divergence  angle  of  incident 
beam. 

0 

3-2 

DELTAE 

Relative  energy  level  of  a 
ray,  with  respect  to  its 
energy  at  the  entrance 
aperture,  below  which  it 
will  be  abandoned. 

10-12 

DEPRT 

Minimum  relative  energy 
level  of  a  ray  for  which  the 
trace  history  will  be  printed 
in  arriving  at  the  exit  aperture. 
Normally  it  is  larger  than 
DELTAE  to  avoid  listing  of  too 
many  rays  with  low  energy. 

10-12 

7-2 

LEMIT 

Number  of  rays  to  start  at  en¬ 
trance  aperture  for  external 
radiation  calculation. 

100 

LREN 

Limit  on  number  of  reflections 
and  transmissions  of  a  ray. 

20 

ISEEDP 

After  every  ISEEDP  rays  are 
finished,  the  value  of  ISEED 
is  printed  out.  These  ISEED 
values  can  be  used  to  continue 
the  computation. 

100 

*  All  angles  are  in  degrees  except  when  otherwise  noted. 
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Default 

Value 


Variable 

NPRINT 

NERR 

NSEC 

RMAP 

KPRINT 

WAVEL 

DEDGE 

DRENT 


Definition 

Limit  on  number  of  rays 
for  which  the  tracing  will 
be  printed  in  arriving  at 
exit  aperture. 


200 


Number  of  error  occurrences  20 
which  will  terminate  program 
execution.  The  accumulated 
energy  distribution  on  the  exit 
aperture  will  be  printed  out. 


Number  of  sections.  0 


Reference 


2-16 


Radius  of  the  energy  distribu-  0  7-3 

tion  maps  at  exit  aperture. 

Signal  for  printing  ray  history  0  7-2 

of  rays  that  didn't  reach  the 
exit  aperture.  Set  any  positive 
integer  and  ray  history  will  be 
printed  out  after  first  KPRINT 
rays.  Set  0  if  no  such  print¬ 
out  is  desired. 

Wave  length  in  microns.  0 

Number  of  wave  lengths  from  100  3-15 

the  edge  of  an  aperture  be¬ 
yond  which  diffraction  will  be 
ignored. 

Width  of  the  outer  annular  ring  0  4-4 

on  the  entrance  aperture  where, 
in  the  case  of  diffraction,  50% 
of  the  rays  will  be  concentrated. 

Set  zero  if  no  such  importance 
sampling  is  desired  at  entrance 
aperture. 
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Variable 

Definition 

Default  Reference 

Value  Page 

ISENT 

Surface  number  of  the 
entrance  aperture. 

0 

ISEXT 

Surface  number  of  the  exit 
aperture. 

0 

AREAEN 

Area  of  the  entrance 
aperture. 

The  area  of  a  circle 
of  radius  PTLCTR 
(7,  ISENT) 

AREAEX 

Area  of  the  exit 
aperture. 

The  area  of  a  circle 
of  radius  PTLCTR 
(7,  ISEXT). 

6.  3  NAMELIST  NAME2 


Diffraction  command,  degrees  of  surface  equations,  and  indices  of 
refraction  are  listed  below. 

Default  Reference 


Variable 


KDIF  (I) 


IDGREE(I) 
RINDEX  (1) 


Definition 

Signal  of  diffraction  at 
apertures.  Set  KDIF  (I) 

=  1  for  I  =  Surface  number 
of  all  the  apertures  where 
diffraction  is  considered. 

Degree  of  equation  for 
surface  I. 

Index  of  refraction  for 
section  I. 


6.4  NAMELIST  NMGEOM 
Basic  surfaces  are  shown  below. 
Variable  Definition 


IP  LANE  (I) 


Surface  number  of  the  Ith 
plane. 


Value 


Page 


3-15 


Given  by 
basic  surfaces. 

1 


2-1 


3-11 


Default 

Value 


Reference 

Page 

2-2 
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Default  Reference 


Variable 

Definition  Value 

Page 

— 

ICONE  (I) 

Surface  number  of  the 

Ith  cone. 

0 

2-2 

ISPHER  (I) 

Surface  number  of  the  Ith 
sphere. 

0 

2-4 

IPARAB  (I) 

Surface  number  of  the  Ith 
paraboloid. 

0 

2-4 

IBAFF  (I) 

Surface  number  of  the  knife 
edges  of  the  Ith  baffle.  Note 
that  there  are  four  surfaces 
for  each  baffle. 

0 

2-4 

PTLCTR  (K,  J) 

Parameters  describing  surface 

J:  The  kind  of  basic  surface 
is  given  by  the  parameters 

IP  LANE,  ICONE,  ISPHER, 
IPARAB  and  IBAFF.  Meaning 
of  PTLCTR  is  fully  described 
in  Section  2. 

0 

2-2 

through 

2-7 

PKNIFE 

Radius  of  the  cross  section  of 
knife  edges. 

0 

2-4, 

2-6, 

2-8 

TBAFF 

Thickness  of  baffles. 

0 

2-4, 

2-8 

2-6, 

DBAFF 

Half  width  of  hyperboloids  on 
plane  baffles. 

0 

2-4, 

2-6 

ABAFFD 

Blade  angle  of  baffle  knife  edges. 

0 

2-4, 

2-8 

2-6, 

COEF  (J.  I) 

Coefficients  of  the  surface 
equation,  Eq{2-1),  for  surface 

0 

2-1 
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6.5  NAMELIST  NAME3 

Special  constraints  are  listed  in  this  subsection. 

Default  Reference 

Variable  Definition  Value  Page _ 

NSTRT  (I)  Number  of  constraints  of  Number  of  2-13 

surface  I.  standard 

constraints. 

KSIDSP  (J,I)  Surface  number  of  the  0  2-13 

surface  to  be  used  as 
the  Jth  constraint  of 
surface  I.  Set  -1  if  no 
surface  is  available  and 
a  new  surface  is  to  be  read 
in  as  the  constraint.  Set  0 

if  the  standard  constraint  is 

used. 

KSTRT  (J,  I)  Type  of  constraint  -  the  Jth  Type  of  the  2-14 

constraint  of  surface  I.  standard  con¬ 
straint  if  co¬ 
efficients  of  the 
surface  has  been 
generated  otherwise  0. 

PTLCSP  (K,  J,  I)  Parameters  of  the  new  surface  0  2-14 

to  be  used  as  the  Jth  constraint 
of  surface  I. 

CSTRT  (K,  J,  I)  Coefficients  of  the  constraint  0  2-10 

surface  equation,  Eq.  (2-2), 
for  the  Jth  constraint  of 
surface  I. 


6.6  NAMELIST  NAME4 

System  Configuration  is  listed  below. 

Variable  Definition 

NSFACE  (I)  Number  of  surfaces  of 

Section  I. 


Default  Reference 

Value  Page _ 

0  2-16 
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Variable 

Default 

Definition  Value 

Reference 

Page 

ISFACE  (J,  I) 

Surface  number  of  the  Jth 
surface  of  Section  I. 

0 

2-16 

ISFTP  (J,  I) 

Surface  type  of  the  Jth 
surface  of  Section  I. 

0 

2-17 

KCOEF  (J,  I) 

Coefficient  indicating 
which  side  of  the  Jth 
surface  Section  I  is  on 

0 

2-18 

IMPSF  (J,  I) 

Parameter  relating  important 
surfaces  to  the  Jth  surface 
of  Section  I. 

0 

2-21 

NIMPS  (I) 

Number  of  important  surfaces 
in  Section  I. 

0 

4-5 

IMPORT  (J,  I) 

Surfaces  number  of  the  Jth 
important  surface  of  Section  I. 

0 

4-5 

NRING  (J,  I) 

Number  of  rings  of  equal  area 
in  important  surface  IMPORT 
(J,  I). 

0 

4-6 

NSEG  (J,I) 

Number  of  angular  segments,  of 
equal  area,  for  each  ring  of  the 
important  surface  IMPORT  (J,  I). 

0 

4-6 

6.  7  NAMELIST  NMCOAT 

Coefficients  of  radiation  properties  are  shown  below. 

Default 

Variable  Definition  Value 

Reference 

Page 

COAT  (J,  I) 

Coefficients  of  radiation 
property  of  surface  I. 

0 

3-7 

I EXP  (I) 

Type  of  near  specular 
reflection  model  of  Surface  I. 

0 

3-9 
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6.  8  NAMELIST  NMEMSN 


Miscellaneous  data  for  internal  thermal  emission  are  listed  in  this 
subsection. 

Default  Reference 

Variable  Definition  Value  Page _ 


NPOLE 


Number  of  poles  where  pilot  0  3-4 

rays  are  generated. 


LEMIT 


Number  of  pilot  rays  leaving  Value  given  in 
each  pole.  namelist  NAME1. 


POLE  (J,  I) 

Co-ordinates  of  the  Ith  pole. 

NSFEM  (I) 

Number  of  emitting  surfaces 
associated  with  the  Ith  pole. 

KSECEM  (J,  I) 

Section  number  of  the  Jth  sur¬ 
face  of  the  Ith  pole. 

KSIDEM  (J.  I) 

Surface  number  of  the  Jth 
surface  of  the  Ith  pole. 

TEMP  (I) 

Temperature  of  the  Ith  surface. 

3-4 

3-5 

3-5 

3-5 

3-5 


STBLTZ 


Stefan- Boltzmann  constant.  3.657E-11  3-3 

watt/ inch  2°K'* 


DELTAE 


Energy  level  of  a  ray  below  Value  given  in 
which  it  will  be  abandoned.  Namelist  NAME1. 


DEPRT  Minimum  energy  level  of  a  Value  given  7-2 

ray  of  which  the  trace  history  in  Namelist 
will  be  printed  in  arriving  at  NAME1. 
the  exit  aperture. 


LRFN  Limit  on  number  of  reflections  Value  given 

and  transmissions  of  a  ray.  in  Namelist 

NAME1. 
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NPRINT 

NERR 


KPRINT 


WVBDL 

WVBDH 
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Definition 

After  every  ISEEDP  pilot 
rays  are  traced,  the  value 
of  ISEED  is  printed  out. 
These  ISEED  values  can  be 
used  to  continue  the  com¬ 
putation 


Default  Reference 

Value  Page _ 

Value  given 
in  Namelist 
NAME1. 


Limit  on  number  of  rays  for  Value  given 
which  the  tracing  will  be  in  Namelist 

printed  in  arriving  at  the  exit  NAME1. 
aperture. 


Number  of  errors  having  Value  given  in 

occurred  when  the  compu-  Namelist  NAME1. 
tation  will  be  terminated.  The 
accumulated  energy  distribution 
on  the  exit  aperture  will  be 
limited  out. 


Signal  for  printing  ray  history  Value  given  7-2 

of  those  that  didn't  reach  the  in  Namelist 

exit  aperture.  Set  any  positive  NAME  1. 

integer  and  ray  history  will  be 

printed  out  after  first  KPRINT 

pilot  rays.  Set  0  if  no  such 

print  out  is  desired. 

Lower  and  upper  limits,  respec-  0 
tively,  of  the  wave  band  in  microns. 

When  WVBDH  is  0,  the  emission 
will  cover  the  whole  spectrum. 
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Section  7 
OUTPUT 


The  output  listing  of  an  execution  is  usually  self-explanatory.  For  con¬ 
venience  to  the  user,  however,  we  list  below  the  complete  contents  of  an 
output  according  to  their  proper  sequence.  The  contents  are  divided  into 
seven  major  groups,  some  of  which  may  not  appear  in  the  result  of  a 
particular  computation.  Except  for  the  error  diagnosis,  we  briefly  describe 
in  this  section  the  significance  of  the  output. 


7.  1  INPUT  DATA  AND  SYSTEM  MODEL 

The  input  data  set  is  first  listed  exactly  as  it  appears  on  the  data  deck. 
After  subsequent  listing  of  the  parameters  in  namelists  NAME1  and 
NAME2,  tables  are  constructed  for  planes,  cones,  spheres,  paraboloids 
and  baffles.  This  is  followed  by  the  numbers  of  constraints  and  a  table 
for  special  constraints. 

The  array  ISFACE  that  defines  the  surfaces  is  given  for  each  section.  The 
corresponding  parameters  ISFTP,  KCOEF,  and  IMPSF  are  listed  in 
similar  tables.  The  important  surfaces,  if  any,  are  then  given  for  each 
section,  together  with  the  numbers  of  rings  and  segments  into  which  each 
important  surface  is  to  be  partitioned. 

In  the  table  of  radiation  properties,  a  set  of  components  is  listed  for  each 
surface.  For  each  of  the  surfaces  with  near  specular  components,  the 
angular  decay  constants  are  given  at  the  end  of  the  table. 

The  ten  coefficients  in  Equation  2-1  are  listed  for  each  of  the  surfaces 
except  toroids.  Immediately  following  these  coefficients,  the  type  of  con¬ 
straint  and  the  coefficients  of  the  constraint  surface  as  expressed  by 
Equation  2-2  are  given  for  each  of  the  constraints.  In  the  present  version 
of  the  program,  Cio  has  the  opposite  sign  that  it  would  have  according  to 
Equation  2  -2. 

The  coefficients  given  for  toroids  are  expressed  in  local  coordinates  and 
are  meaningless  to  the  user  since  the  equation  has  not  been  presented  in 
the  manual.  There  is  no  input  data  directly  involved,  however,  and  the 
user  can  disregard  these  coefficients. 
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7.  2  EXTERNAL  RADIATION  RAY  TRACING 

When  a  ray  reaches  the  exit  aperture  and  the  relative  energy,  with  respect 
to  its  energy  when  it  started  at  the  entrance  aperture,  is  greater  than 
DEPRT,  the  trace  of  the  ray  will  be  printed  out  as  in  the  following: 


KEMIT  KPOLE  En 


where  KEMIT  is  the  ray  number,  KPOLE  the  pole  number,  and  En  is  the 
relative  energy  of  ray  when  it  read  es  the  exit  aperture.  Since  poles  are 
not  used  in  external  radiation,  KPOLE  will  he  zero  in  this  part.  (Xj,  Yj, 
Zj)  is  the  initial  location  of  the  ray  at  the  entrance  aperture,  (X2,  Y2,  Z2) 
is  the  first  intercept,  and  so  forth.  (Xn,  Yn,  Zn)  is  the  intercept  at  the 
exit  aperture.  E^  is  the  relative  energy  of  the  ray  at  the  initial  point, 

E2  the  relative  energy  at  the  first  intercept,  and  so  forth. 

When  KPRINT  is  a  positive  integer,  traces  are  printed  out  for  all  rays 
with  KEMIT  >  KPRINT  regardless  of  whether  they  reach  the  exit  aperture 
or  not. 

After  each  group  of  ISEEDP  rays  has  been  initiated  and  finished,  KEMIT 
and  ISEED  are  listed.  By  using  one  of  these  ISEED  values  as  the  initial 
ISEED  for  the  next  computation,  the  user  can  continue  calculation  at  some 
intermediate  point  or  at  the  end.  This  is  useful  in  troubleshooting  in  the 
moddle  of  a  calculation. 

Error  messages  appearing  in  this  part  will  be  discussed  in  Section  8. 
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7.  3  RAY  COUNT  AND  ENERGY  DISTRIBUTIONS  OF  EXTERNAL  RADIATION 

a.  Spatial  Ray  Count  Distribution 

A  square  of  sides  2*  RMAP  is  divided  into  a  grid  of  30  by  30  square  elements. 
The  horizontal  distance  corresponds  to  the  local  x-axis  at  the  exit  aperture 
plane  while  the  vertical  distance  the  local  y-axis.  The  number  at  each 
element  indicates  the  number  of  rays  striking  the  element. 

b.  Spatial  Energy  Distribution 

Energy  distribution  is  printed  out  in  the  same  way  as  the  ray  count.  Due 
to  limitation  in  space,  it  is  divided  into  three  portions.  The  first  lists 
the  first  10  columns,  the  second  lists  the  middle  10  columns  and  the  rest 
appears  on  the  third.  The  values  are  based  on  a  unit  energy  input  at  the 
entrance  aperture. 

c.  Radial  Ray  Count  Distribution 

The  circle  of  radius  RMAP  at  the  exit  aperture  is  divided  into  100  rings  of 
equal  width.  Listed  here  are  numbers  of  rays  in  each  of  the  rings. 

d.  Radial  Energy  Density  Distribution 

Similar  to  "c",  the  energy  density  distribution,  energy  per  unit  area  per 
unit  energy  input,  is  listed  for  the  100  rings.  The  energy  density  distri¬ 
bution  is  then  plotted  against  the  radial  distance  in  the  figure  that  follows. 

e.  Spatial  and  Angular  Ray  Count  Distribution 

The  square  of  side  2*  RMAP  on  the  exit  aperture  is  divided  into  a  grid  of 
3  by  3  area  elements  as  follows: 
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where  x  and  y  are  the  local  coordinate  axes.  The  numbers  1  through  9 
refer  to  the  elements.  The  rays  are  grouped  by  the  elements  they  hit  in 
arriving  at  the  aperture  plane.  In  addition,  they  are  grouped  by  the  angles 
between  the  ray  and  the  exit  aperture  normal. 

f.  Spatial  and  Angular  Energy  Distribution 

Similar  to  "e",  energy  distributions  are  listed  according  to  the  elements 
and  the  angles  to  the  normal  when  they  reach  the  exit  aperture. 


7.  4  SUMMARIZED  RESULTS  OF  EXTERNAL  RADIATION 

a.  Number  of  Rays  Entering 

Number  of  rays  actually  started  at  entrance  aperture.  This  will  be  the 

same  as  LEMIT  except  when  computation  is  terminated  due  to  excessive  errors. 

b.  Rays  and  Energy  Through  Entrance  Aperture 

Number  of  rays  and  the  total  relative  energy  that  leaves  the  system  via 
the  entrance  aperture. 

c.  Rays  Through  Exit  Aperture 

Number  of  rays  hitting  the  exit  aperture.  There  are  three  numbers.  The 
first  covers  all  the  rays  in  the  first  half,  that  is,  REMIT  <  LEIvIIT/2,  the 
second  cover  those  in  the  second  half,  and  the  last  is  the  total  number. 

The  number  is  divided  into  two  parts  for  an  indication  of  the  accuracy. 

Since  a  statistical  analysis  of  the  error  is  also  given,  the  user  could  dis¬ 
regard  this  division. 

d.  Total  Energy  Through  Exit  Aperture 

Energy  reaching  the  exit  aperture  per  unit  energy  entering  the  system. 

The  first  number  is  from  the  first  half,  the  second  from  the  second  half, 
and  the  last  from  the  whole.  Again,  the  user  could  disregard  the  first  two 
numbers. 
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e.  Attenuation 

Ratio  of  energy  intensity  of  the  incident  beam  to  the  energy  density  at  the 
exit  aperture.  The  latter  includes  energy  bundles  in  all  directions.  The 
first  number  is  the  attenuation  obtained  from  the  first  half,  the  second 
from  the  second  half,  and  the  last  from  the  whole.  Again  the  user 
disregard  the  first  two  numbers. 

f.  Off-Axis  Rejection 

Ratio  of  energy  entering  the  system  to  that  reaching  the  exit  aperture. 

g.  Standard  Deviation /n/ti  •  Average 

Relative  error  of  total  energy  hitting  the  exit  aperture  with  90  percent 
confidence. 

h.  Total  Energy  Reaching  Target 

Range  of  the  total  relative  energy  in  "d"  with  90  percent  confidence. 

i.  Attenuation 

Range  of  attenuation  in  "e"  with  90  percent  confidence. 

j.  Maximum  Number  of  Rays  Stored 

The  maximum  number  of  rays  stored  due  to  splitting  at  any  time  of  the 
entire  external  radiation  computation.  The  limit  is  100,  although  it  can 
readily  be  increased  if  necessary.  If  at  any  time  the  number  of  rays 
stored  exceeds  the  limit,  the  exceeding  rays  are  abandoned. 


7.  5  INTERNAL  EMISSION  RAY  TRACING 

After  a  listing  of  the  parameters  in  namelist  NMEMSN,  ray  tracing  is  given 
for-  the  internally  emitted  rays  in  the  same  manner  as  that  given  for  the 
external  rays  described  in  Subsection  7.2.  The  energy  is  in  terms  of  the 
basic  power  unit,  watts. 

Error  messages,  if  any,  will  appear  in  this  part  of  the  output.  The  dis¬ 
cussion  will  be  given  in  Section  8. 
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7.  6  RAY  COUNT  AND  ENERGY  DISTRIBUTIONS  OF  INTERNAL  EMISSION 

The  ray  count  and  energy  distributions  are  given  for  internal  emission  in 
the  same  way  as  those  given  for  external  radiation  (see  Subsection  7.  3). 
Again,  the  energy  is  in  terms  of  the  basic  power  unit. 


7.  7  SUMMARIZED  RESULTS  OF  INTERNAL  EMISSION 

a.  Number  of  Pilot  Rays 

Number  of  pilot  rays  generated  at  each  pole.  It  is  given  by  LEMIT  except 
when  computation  is  terminated  due  to  excessive  errors. 

b.  Rays  and  Energy  Through  Entrance  Aperture 

Number  of  rays  and  the  total  emitted  energy  that  leave  the  system  via  the 
entrance  aperture. 

c.  Rays  Through  Exit  Aperture 
See  "c"  of  Subsection  7.  4. 

d.  Total  Energy  Through  Exit  Aperture 
Emitted  energy  that  reaches  the  exit  aperture. 

e.  Standard  Deviation/Vn  •  Average 

Relative  error  of  total  emitted  energy  hitting  the  exit  aperture  with  90 
percent  confidence. 

f.  Total  Energy  Reaching  Target 

Range  of  the  total  emitted  energy  in  "d"  with  90  percent  confidence. 

g.  Maximum  Number  of  Rays  Stored 

The  maximum  number  of  rays  stored  due  to  splitting  at  any  time  of  the 
entire  internal  emission  computation.  This  limit  is  100,  although  it  can 
readily  be  increased  if  necessary.  If  at  any  time  the  number  of  rays  stored 
exceeds  the  limit,  the  exceeding  rays  are  abandoned. 
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Section  8 
DIAGNOSTICS 

Aside  from  minor  errors  in  a  Namelist  data,  for  example,  a  mispunch 
of  a  character,  most  of  the  errors  occur  in  the  middle  of  the  ray  tracing. 
This  section  presents  such  error  messages  and  the  procedure  to  follow 
in  tracking  down  the  sources  of  errors.  Each  error  message  is  accom¬ 
panied  by  a  set  of  informative  parameters.  Definitions  are  first  given 
below  for  the  frequently  used  parameters. 

KEMIT  Hay  number,  from  1  to  LEMIT 
IPOLE  Pole  number,  from  1  to  NPOLE 

K  STY  PE  Type  of  surface,  given  by  ISFTP  of  the  surface  that 

intercepts  the  ray  immediately  before  the  error.  When 
no  interception  of  a  ray  appears  after  testing  with  all 
the  surfaces  in  the  section,  KSTYPE  is  set  equal  to  6. 

KSID  Surface  number  of  the  surface  of  latest  interception. 

E  Energy  of  the  ray 

External  Radiation  Per  unit  energy  input  at  entrance 

aperture 

Internal  Emission  In  absolute  power  units 

RAY  Location  (X,  Y,  Z)  and  direction  cosines  (a,  p,  y)  of  the 

ray 

KSEC  Section  number 

ISF  Sequence  number  of  the  important  surface  in  the  section. 

Note  that  there  can  1  j  up  to  three  important  surfaces  in 
each  section. 

SURNOR  Direction  cosines  of  the  inward  unit  surface  normal  at 
the  point  of  interception. 

a.  Error  in  Emission  -  KSID,  IPOLE.  Tj,  T^,  RAY 

Calculation  of  the  interception  between  the  pilot  ray,  location  and  direction 
given  by  RAY,  and  surface  KSID  results  in  the  equation: 

T0  +Tj  S=  0 
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where  S  is  the  distance  from  the  pole  to  the  point  of  interception,  and  Tg 
and  T j  are  parameters  obtained  from  substitution.  The  fact  that  the 
magnitude  of  T  j  approaches  zero,  less  than  10"^,  renders  the  equation 
indeterminate.  The  user  should  examine  whether  the  surface  is  parallel 
to  the  pilot  ray.  The  problem  can  then  be  resolved  by  relocating  the  pole. 

b.  No  Intercept  -  KEMIT,  K  STY  PE,  KSEC,  KS1D,  E,  RAY 

Each  section  is  a  complete  enclosure.  When  no  interception  can  be  found 
for  a  ray,  either  some  surface  is  missing  or  errors  exist  in  the  surface 
or  the  constraint  parameters. 

If  a  quick  check  fails  to  reveal  the  source  of  error,  the  user  is  advised  to 
examine  the  ray  location  and  direction  in  terms  of  the  section  configuration. 
If  the  ray  is  correctly  situated  in  the  section,  he  should  then  observe  the 
section  and  find  out  the  surface  that  should  have  intercepted  the  given  ray. 
The  coefficients  of  the  surface  and  its  constraint  equation  (see  Subsection 
7.  1)  will  pinpoint  the  error  source. 

The  user  should  also  check  the  input  values  of  NSFACE,  ISFACE,  ISFTP 
and  KCOEF  (see  Subsection  6.6). 

When  the  above  error  occurs,  the  ray  is  abandoned.  However,  when  the 
number  of  errors  exceeds  NERR,  the  computation  is  terminated. 

c.  Error  in  GEOM  -  KSTYPE 

The  value  of  KSTYPE  is  limited  1  through  5.  Examine  ISFTP  portion  of 
the  input  data  for  possible  error. 

d.  Error  in  GEOM,  THIRD  DEGREE,  polynomial  is  not  available  -  KSID, 
IDGREE. 

The  value  of  IDGREE  is  limited  to  1,  2  and  4.  Examine  IDGREE  portion 
of  the  input  data  for  possible  error. 

e.  Error  in  NORM.  The  Third  Degree  Polynomial  not  available. 

Examine  IDGREE  input  for  possible  error. 
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f.  Secondary  Ray  Trapped  by  Important  Surface  in  REFLN  -  KEMIT, 

I  POLE,  KSEC,  KSID,  ISF,  RAY,  SURNOR 

After  an  important  ray  is  generated  and  stored,  the  direction  of  the  random 
ray  which  carries  the  remaining  energy  is  obtained  according  to  the  parti¬ 
cular  model  of  energy  density  distribution.  According  o  the  theory  how¬ 
ever,  the  random  ray  should  not  hit  the  important  surface.  If  the  ray  so 
chosen  happens  to  hit  the  important  surface,  a  new  direction  is  selected. 
This  process  is  repeated  until  a  ray  is  obtained  that  misses  the  important 
surface . 

To  avoid  unnecessary  waste  of  computation  time,  a  ray  is  abandoned  after 
it  fails  to  miss  the  important  surface  seven  times.  The  diagnostic  is 
printed  out  to  assist  in  deciding  whether  the  important  surface  set  should 
be  modified. 

g.  Error  in  REFRN  -  KSID,  KSEC,  KSTYPE,  IRTEST,  COSPSN 

The  angle  between  the  ray  and  the  normal  should  be  between  0  and  n/2. 

If  it  is  outside  this  range,  the  error  message  is  printed  out  and  the  angle 
is  defaulted  to  0  or  ?r/2  according  to  which  side  the  value  falls  on. 

The  parameter  IRTEST  is  1  in  case  of  total  internal  reflection. 

COSPSN  is  the  cosine  of  the  angle  between  the  ray  and  the  inward  surface 
normal. 

h  Error  in  ANGSSD  -  PSI  (angle  between  the  specular  beam  and  the 

normal),  ALPHA  (the  angle  between  the  reflected  ray  and  the  specular 
reflected  beam),  BETA  (the  azimutnal  angle  of  near  specular  scattering), 

ALPHAM  (maximum  angle  of  ALPHA)  and  PHI  (the  angle  between  the 
reflected  ray  and  the  normal). 

Angle  PHI  should  be  between  0  and  n/2.  When  it  is  outside  this  range,  the 
message  is  printed  out  to  warn  of  the  error.  The  ray  will  be  abandoned  in 
the  subsequent  calculation  in  subroutine  GEOM  when  no  interception  can  be 
obtained. 

i.  Error  in  SCATTR,  no  solution  found  for  X  -  last  value  equals  X  and 
difference  DR  equals  DR, 

Solution  of  the  equation 

X  *  f(a). 
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which  is  obtained  from  the  near  specular  terms  of  Equation  3-1,  by  iteration, 
does  not  converge,  X  is  the  probability  function  and  a  is  the  off  specular  angle. 

j.  Error  in  DIFRN  -  KSEC,  RAY  and  DELA  (the  distance  from  the  inter¬ 
secting  point  to  the  aperture  edge) 

An  interception  of  a  ray  by  an  aperture  is  determined  before  computation 
enters  subroutine  DIFRN.  Thus  the  shortest  distance  from  the  intersecting 
point  to  the  aperture  should  always  be  greater  than  zero.  When  it  is  zero 
or  negative,  ihe  message  is  printed  out  to  warn  of  possible  error  in  the 
input  data  of  the  surfaces  involved.  The  ray  will  be  traced  without  being 
diffracted  at  this  point. 

k.  Secondary  Ray  Trapped  by  Important  Surface  in  DI^RN  -  KEMIT, 

I  POLE,  KSEC,  KSID,  ISF,  RAY,  SURNOR,  E. 

See  "f". 

l.  Negative  Energy  Ray  -  E,  ESAVE  (part  of  E  that  belongs  to  the 
important  ray) 

FACl  (consine  of  the  angle  between  the  important  ray  and  the  unit  normal 
of  the  important  surface),  FAC2  (factor  of  the  energy  going  in  the  direction 
of  the  important  ray),  DENOM2  (  quare  of  the  length  of  the  important  ray), 
KEMIT.  KSTYPE,  KSEC,  KSID,  ISF,  first  half  of  RAY  (initial  point  of 
the  important  ray),  POINT  (end  point),  OQ  (vector  of  the  important  ray) 
and  SURNOR. 

When  ESAVE  is  negative  in  subroutine  EXPECT,  the  program  by-passes 
splitting  the  ray.  When  FACl  is  positive,  it  is  in  most  cases  due  to  im¬ 
proper  orientation  of  the  important  surface.  The  ray  has  come  from  the 
back  side  of  the  important  surface.  Sometimes,  it  may  be  necessary  to 
arrange  the  important  surfaces  such  that  this  does  occur.  In  this  case, 
the  computation  result  is  not  affected. 

m.  Error  in  EXPECT  -  PHI  (angle  between  the  important  ray  and  the 
surface  normal),  CBETA  (cosine  of  the  azimuthal  angle  of  scattering), 
ALPHAM  (maximum  angle  from  the  specular  beam  to  the  surface 
tangent  for  a  particular  BETA) 

The  important  ray  of  near  specular  scattering  is  transmitted  into  the  surface 
instead  of  being  reflected.  Subsequent  calculation  will  show  no  interception 
and  the  ray  will  then  be  abandoned.  See  Figure  3-3, 
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n.  Excessive  Energy  Given  to  Important  Rays  in  EXPEC1  -  KEMIT, 

I  POLE,  KSEC,  KSID,  ISF,  E,  RAY  (X,  Y.  Z) 

The  important  surface  is  partitioned  and  NRING  x  NSEG  important  rays 
are  generated.  Ideally,  the  total  amount  of  energy  associated  with  these 
important  rays  is  less  than  the  total  energy  leaving  the  surface.  However, 
since  the  energy  of  each  important  ray  is  approximated  by  using  a  randomly 
selected  point  in  each  partition,  their  sum  might  exceed  the  actual 
total  energy  leaving  the  surface. 

When  the  energy  is  exhausted  before  NRING  x  NSEG  important  rays  are 
generated,  the  remaining  important  rays  and  the  complementary  ray  are 
abandoned.  The  user  should  examine  the  important  surface.  As  we  see 
from  Section  4,  the  situation  can  be  improved  by  reducing  the  size  of  the 
partitions,  that  is,  increasing  NRING  and/or  NSEG. 

o.  Negative  Energy  Ray  in  DIFEXP  -  FAC1  (cosine  of  the  angle  between 
the  important  ray  and  the  unit  normal  of  the  important  surface),  DENOM 2 

(square  of  the  length  of  the  important  ray),  KEMIT,  KSTYPE,  KSEC, 
ISF,  RAY,  OQ  (vector  of  the  important  ray)  and  POINT  (end  point) 

When  ESAVE  is  negative  in  subroutine  DIFEXP,  program  by-passes  splitting 
of  ray.  Very  likely  the  error  is  due  to  incorrect  orientation  of  the  important 
surface.  See  "1". 

p.  Negative  Energy  Ray  in  DIFEXP-E,  ESAVE,  FAC1,  DENOM 2  (see  V), 

A  (^  [(^4)2  +  (fli)2]  ),  TTHES  (tan  0*)t  TTHIS  (tan  0*),  RRB  (area  of 

^  o  12 

1  £ 

the  important  surface  hr),  KEMIT,  KSTYPE,  KSEC,  ISF,  RAY,  OQ, 
POINT  (see  "o") 

When  ESAVE  is  negative  in  subroutine  DIFEXP,  program  by-passes  splitting 
of  ray.  The  users  should  examine  the  parameters  WAVEL,  DEDGE  and 
DRENT  and  the  aperture  surface  through  which  the  diffraction  occurs. 

q.  Excessive  Energy  Given  to  Important  Rays  in  DIFEXP  -  KEMIT, 

IPOLE ,  KSEC,  KSID,  ISF,  E,  RAY  (X,  Y,  Z) 

Similar  to  that  of  EXPECT  given  in  "n". 
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r.  Error  in  TESTSS  -  KEMIT,  KSEC,  KSID,  ISF,  IMPS  (surface  number 
of  the  important  surface),  Tj,  Tq  (the  two  coefficients  in  the  equation 
for  the  distance  of  the  ray  to  the  important  surface)  and  RAY. 

The  splitting  point  is  too  close  to  the  extended  important  surfaces.  The 
testing  of  the  secondary  ray  of  being  intercepted  by  the  important  surface 
is  ignored.  The  user  should  check  to  see  if  the  important  surfaces  need  to 
be  modified. 

s.  IRAY  Exceeds  NRAY  -  KEMIT,  KSTYPE,  KSEC,  KSID,  E,  and  RAY. 

There  is  a  limit  of  100  rays  that  can  be  stored  at  any  time.  To  raise  this 
limit,  the  statement  NRAY  =  100  in  main  program  GUERAP  is  changed  to 
NRAY  =  any  desired  value. 

t.  Error  in  DISTB  -  K  (group  number  with  respect  to  angle),  L  (regions), 
KEMIT,  KSTYPE,  KSEC,  KSID,  RAY,  and  x,  y,  z  (location  in  terms 
of  the  local  coordinates). 

When  RMAP  is  smaller  than  the  exit  aperture  radius,  some  of  the  rays 
that  reach  the  exit  aperture  may  miss  the  ray  count  and  energy  distribution 
maps.  The  rays  are  simply  ignored. 

u.  Error  in  ACOS  -  X  (argument)  and  KEMIT. 

Absolute  value  of  the  argument  of  an  arc -cosine  is  greater  than  unity. 
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Section  9 
EXAMPLES 


A  series  of  examples  are  presented  in  this  section  to  demonstrate  the 
application  of  the  GUERAP  program  in  a  wide  variety  of  problems.  The 
user  is  advised  to  follow  through  these  examples  before  attempting  the 
task  of  simulating  complicated  systems. 

The  first  example  is  a  brute  force  Monte  Carlo  solution  of  a  simple  dif¬ 
fraction  problem.  The  radiant  intensity  distribution  as  a  result  of  diffrac¬ 
tion  effects  is  obtained  on  a  target  screen.  Comparison  of  the  result  with 
available  solutions  is  presented  in  more  detail  in  Diffraction  Modeling 
(Reference  2). 

When  the  radiant  intensity  is  sought  in  only  a  small  portion  of  the  screen, 
the  importance  sampling  technique  is  applied  with  great  success.  This  is 
presented  in  the  second  example. 

Examples  3  and  4  deal  with  internal  thermal  emission  of  simple  cavities. 
The  result  of  the  third  example  is  compared  with  an  available  independent 
solution. 

In  the  last  example,  we  compute  the  attenuation  of  a  monochromatic  colli¬ 
mated  beam  through  a  typical  optical  system.  The  example  makes  use  of 
various  features  available  in  the  program  such  as  knife-edge  baffles, 
lenses  and  mirrors.  It  also  demonstrates  the  use  of  importance  sampling 
in  an  optical  system  with  high  attenuation. 


9.  1  EXAMPLE  1 

For  the  first  example,  we  consider  the  simple  system  shown  in  Figure  9-1. 
A  laser  beam  of  10.6  p  wave  length  enters  a  circular  cylinder  of  3  inches 
diameter  through  a  small  circular  aperture  of  0. 1  inch  diameter.  The 
aperture  has  a  perfect  knife-edge  so  there  will  be  no  reflection  off  the 
aperture  edge.  The  energy  is  collected  on  a  screen  100  inches  behind  the 
entrance. 

In  the  absence  of  diffraction,  the  energy  will  be  uniformly  distributed  over 
a  circle  of  the  same  size  as  that  of  the  entrance  aperture.  When  diffraction 
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FIGURE  9-1.  SCHEMATIC  OF  DIFFRACTION  TEST 

effects  are  included,  however,  the  energy  distribution  on  the  screen  should 
show  an  approximation  of  a  Fraunhofer  diffraction  pattern  (Reference  4,  5), 
This  example  was  originally  used  to  check  the  diffraction  model  of  the 
program. 

There  are  four  surfaces,  three  planes  and  one  cone,  in  the  system  when 
importance  sampling  is  not  utilized.  The  boundaries  of  all  the  surfaces  can 
be  completely  described  by  standard  constraints  so  no  special  constraint  is 
needed.  The  complete  set  of  input  data  is  given  in  Table  9-1. 

The  title  card  is  terminated  by  ENDTITLE  starting  at  column  1.  This  is 
followed  by  ISEED  =  1  in  020  format. 

The  laser  beam  enters  the  system  parallel  to  the  Z-axis  and  is  assumed 
completely  collimated,  so  SUNAGD,  SUNAZD  and  SUNDVD  are  all  zero,  as 
given  by  default. 

Energy  levels  of  a  ray  for  termination  and  listing  of  ray  tracing,  DELTAE 
and  DEPRT,  are  both  at  10'  ^  by  default.  This  is  immaterial  because  the 
rays  will  either  reach  the  screen  directly  or  be  absorbed  completely  by  the 
black  body  inner  surfaces. 
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TABLE  9-1.  INPUT  DATA  OF  EXAMPLE  1 


DIFFRACTION  TESTING 
NO  IMPORTANT  SAMPLING 

ENEKoY  DENSITY  DISTRIBUTION  SOUGHT  FOR  ENTIRE  TARGET  SCREEN 

ENUTITLE 

1 

SNAME1  LEMIT  =  20000»LRFN=1*  TSEEDP=100»NSEC=i *KMAP=1 .5*rtAVEL  =  10.6% 
UWENTsO»ISENT=J«ISEXT=4«KPR INT  =  0  *DEDGF=<?004 
4NAME2  KOIF ( 3) =1 i 
4NMGEOM  IPLANE ( 1 >  =2»  3«4« 

PTLCTRI 1*2) =0,0.1  00 ♦ 0 . 0 . -1  , I . 5* . 05  ♦ 

PTLCTR( I«3)=0«0t  100«0*0«-1  ♦.05*0» 

PTLCTW(l»4)=OtO»0*0*0«l»l«StO* 

ICONE ( l ) =1 « 

PTLCTR(l,I)s0«0tl00*0t0t0«1.5*l.S4 
4NAMF3  4 

$NAME4  NSFACE(1)=4,  I SFACF < 1 ♦ 1 > *1 *2» 3*4,  KCOEF  ( 1 « 1  >*<)•  1  •  1  •  1  • 

ISFTP(1»1)=S»5«1«2S 

SNMCOAT  COAT ( 1 « 1 )  =  1  *  C0AT(1*2)-14 

4NMEMSN  LEM  I T  =  04 
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A  small  number  of  rays  are  needed  for  a  simple  system  when  we  are 
looking  for  the  overall  attenuation.  However,  a  large  number  is  required 
if  we  are  interested  in  an  accurate  energy  distribution  over  the  screen. 

Thus  20000  rays  are  generated.  The  limit  on  the  number  of  reflection 
LFRN  is  set  at  1  because  we  axe  collecting  only  the  energy  that  reaches 
the  screen  directly  from  the  entrance  aperture.  [SEED  is  to  be  printed 
once  every  100  rays. 

By  not  specifying  a  value  for  either  NPRINT  or  NERR,  the  number  of  rays 
to  be  printed  out  is  200  and  the  maximum  number  of  errors  permitted  is 
20  by  default.  There  is  only  one  section  in  the  system,  thus  NSEC  ■  I. 

The  energy  distribution  is  wanted  over  the  entire  target  screen,  thus  we 
set  RMAP  at  1.5  inches.  K PRINT  is  zero  for  the  first  trial  when  no  addi¬ 
tional  diagnostics  are  needed.  The  wave  length  of  the  incident  beam  is 
10. Gp.  The  example  is  mainly  dealing  with  diffraction  at  the  entrance 
aperture,  so  it  is  essential  to  include  the  entire  aperture  area  in  our 
consideration  of  diffraction.  This  is  achieved  by  assigning  200  to  the  value 
of  DEDGE,  the  number  of  wave  lengths  necessary  to  cover  the  entire  dif¬ 
fraction  zone.  No  importance  sampling  is  needed  at  the  entrance  aperture 
so  DRENT  is  zero.  Surfaces  3  and  4  are,  respectively,  the  entrance  and 
the  exit  apertures. 

Diffraction  is  sought  for  the  entrance  aperture  with  surface  number  3,  so 
we  set  KDIF(3)  =  1. 

The  surface  number  of  the  planes  are  2,  3  and  4,  and  that  of  the  cone  is  I. 
Referring  to  Figure  2-1,  we  see  that  surface  2  can  be  specified  by  the  center 
of  aperture  (0,  0,  100)  and  the  unit  normal  (0,  0,  -1),  which  is  to  be  preferred 
pointing  into  the  system,  although  not  mandatory.  The  surface  is  constrained 
by  two  circular  cylinders  around  the  prescribed  center  point  and  unit  normal. 
The  outer  constraint  has  a  radius  of  J.5  while  the  inner  one  0.05.  These  are 
the  sequence  of  numbers  specified  by  the  array  PTLCTR  for  the  surface. 

Surfaces  3  and  4  are  similarly  specified  except  that  there  are  no  inner  con¬ 
straints.  The  cone  surface  is  described  by  the  two  center  points  (0,  0,  100) 
and  (0,  0,  0),  and  the  radii  on  both  ends,  1.5  and  1.5;  sec  Figure  2-2. 

These  are  given  in  the  array  PTLCTR  for  the  surface. 

The  index  of  refraction  is  1  for  air.  Referring  again  to  Figure  9-1,  we 
see  that  there  are  four  surfaces  in  the  section,  namely  1,  ?,  3  and  4.  The 
parameters  PTLCTR  show  that  all  the  unit  normals  for  the  planes  point 
into  the  region,  so  KCOEF  is  1  for  all  the  planes.  It  is  0  for  the  cylinder 
because  the  region,  that  is,  the  section,  is  inside  the  cylinder. 
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The  first  two  surfaces  in  the  sec  ion  are  non-transmitting,  the  third  one 
is  an  aperture  toward  the  entrance  and  the  last  surface  is  an  aperture  toward 
the  exit.  These  are  so  specified  by  array  ISFTP. 

Surfaces  1  and  2  are  black  bodies  with  100  percent  absorptance,  giving 
COAT(l.  1)  =  COAT(l ,  2)  =  1. 

The  resultant  radiant  energy  density  distribution  on  the  target  screen  is 
plotted  against  the  distance  from  the  center  in  Figure  9-2.  For  comparison, 
the  average  densities  at  three  small  radial  bands  are  listed  in  the  second 
row  of  Table  9-3, 


9.2  EXAMPLE  2 

In  this  example,  we  shall  demonstrate  the  use  of  importance  sampling  and 
the  procedure  of  repeating  calculations  on  the  same  system  with  minor 
modifications. 

In  Table  9-3,  the  energy  densities  at  three  narrow  radius  bands  are  listed 
from  the  previous  example.  We  will  show  here  how  to  obtain  these  densities 
with  importance  sampling. 

Table  9-2  lists  the  complete  input  data  for  this  example.  There  are  three 
parts  in  the  data  set,  each  for  one  calculation.  The  first  part  is  a  variation 
of  the  input  data  of  Example  1  given  in  Table  9-1. 

In  order  to  calculate  the  energy  density  in  the  center  circle  and  rings  given 
in  Table  9-3,  we  generate  these  surfaces  as  surfaces  5,  6  and  7.  Since  all 
these  surfaces  are  part  of  surface  4,  the  first  six  parameters  are  identical 
to  those  of  surface  4. 

Surface  5  is  the  center  circle  with  0.015  inch  outer  radius  and  no  inner 
radius.  Surface  6  is  the  ring  with  0.375  inch  outer  radius  and  0.36  inch 
inner  radius.  Similarly,  surface  7  is  the  ring  with  0.75  inch  outer  radius 
and  0.735  inch  inner  radius. 

The  important  surfaces  and  their  assignment  are  given  in  NAME4.  There 
is  one  important  surface,  surface  5,  in  Section  1.  The  surface  is  not  parti¬ 
tioned,  so  NRING  and  NSEG  are  both  unity. 

The  splitting  of  rays  toward  the  important  surface  given  above  is  called  for 
when  rays  pass  the  entrance  aperture,  surface  3.  The  corresponding  number 
of  IMPSF  is  given  a  value  unity,  assigning  important  surface  number  1  to  the 
surface,  and  it  is  zero  for  the  rest  of  the  surfaces. 
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TABLE  9-2.  INPUT  DATA  OF  EXAM  PLE  2 


DIFFRACTION  TESTING 

CENTER  CIRCLE  IMPORTANT  SAMPLING 

ENOTITLE 

1 

SNAMEI  lEM I T  =  500 ♦  L«FN=1.  ISEEDP  =  100»  NPRINT=100*  NSEC=1*  RMAP=1.5. 
*AVEL=10.6t  ORENT  =U «  IStNT  =  3.  ISEXT=4,  KPRlNUOt  DEDGE*200S 
SNAME ?  K0IF(3)  =  U 
4NMGFOM  IPLANt  <  1  )=<?♦ 

PTLCTW(l,2)s0.0*l00»0.0t-ltl.S».0St 
PTLCTRI  1.3)  =£>«  0.1 00*0*0  Ot 

PTLCTR(I*4>*0»0«0*0*0«I*l»e>»0« 

PTLCTWU.5)s0*0»0*0t0*i».0lb*0* 

PTLCTRI 1 .6)*0»0.0»0«0» 1« .375. .36t 
PTLCTw(l,7)=0.0.0»0.0.1..7S..735» 

I  CO  NE  ( 1  >  s  1  * 

PTLCTR(ltl)*0.0»100»0*0*0#l.S*l.S$ 

SNAME3  » 

SNAME*  NSFACtU  )=4,  ISF  ACf  (1*1)  =  1»2*‘>»4«  KCOEF  <1  «  I)  «0 . 1  •  1  *  1  • 

IMPSF  ( 1  •  |  )  sQ  1 0 1 1  *0  •  N I  mPS  ())£)«  IMPORT  ( 1  *  1 )  *St  N«ING(  1 1 1  >  *1  t  NSEG  <  I  .  1  >  *=  1  ♦ 
ISFTP(ltl»sS.S.1.2* 

SNMCQAT  COAT  (  1 » 1 )  =  1  *  C0AT(1*2)=1S 

tNMEMSN  LF.MIUOS 

OIFFRACTION  TESTING 

SECOND  RING  IMPORTANT  SAMPLING 

ENOTITLE 

1 

SNAMEI  LEMIT*500i 
JNAME2  S 
SNMGEOM  % 

SNAME3  S 

INAME4  IMPORT ( 1 .lIsGt 
4NMCOAT  S 
4NMEMSN  LEMITsO* 

DIFFRACTION  testing 

THIRD  RING  IMPORTANT  SAMPLING 

ENOTITLE 

1 

SNAMEI  LEMITsbOO* 

SNAME?  i 
4NMGEOM  % 

SNAME  I  S 

SNAME4  ImPOR  T ( I « 1 ) *7S 
SNMCOAT  S 
SNMEMSN  LEMITsOS 


1 
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TABLE  9-  3.  COMPARISON  OF  ENERGY 
DENSITIES  AT  SELECTED  RING 


Energy  Densities 

CDC-6600 

E 

(watts/in^/watt) 

Execution  Time  I 

(inches) 

0-0.015 

0.360-0.375 

0.735-0.750 

(sec) 

EXAMPLE  1 

No  Importance 

Sampling 

1.981 

0.507 

0.083 

63.2 

EXAMPLE  2 

I*-  portance 

Sampling 

2.033 

11.4 

Imoortance 

Sampling 

0.543 

9.1 

Importance 

Sampling 

0.082 

9.4 
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Notice  that  the  number  of  rays  given  for  the  external  radiation  in  NAME1 
is  500  instead  of  20,000  in  Example  1.  In  Example  1,  the  rays  are  allowed 
to  select  their  directions  randomly  according  to  the  probability  density 
functions  and  the  number  of  rays  to  strike  a  particular  radial  band  such  as 
the  center  circle  is  in  the  order  of  the  total  number  of  rays  divided  by  the 
number  of  bands.  In  order  to  obtain  a  significant  number  of  rays  for  each 
band,  it  was  necessary  to  generate  an  extremely  large  number  of  rays  at 
the  entrance  aperture. 

In  the  first  part  of  Example  2,  each  of  the  rays  generated  is  split  into  two 
rays,  one  of  which,  the  important  ray,  is  aimed  at  the  center  circle.  Thus, 
a  substantially  smaller  number  of  rays  are  required  to  get  good  statistical 

accuracy. 

After  the  calculation  of  energy  density  at  the  center  circle,  we  now  proceed 
to  calculate  the  density  at  the  two  rings,  namely  surfaces  6  and  7.  In  the 
second  and  the  last  parts  of  the  data  set,  we  direct  the  calculations  with 
the  given  modifications.  The  second  part  is  for  surface  6,  while  the  last 
is  for  surface  7.  Note  that  it  is  necessary  to  respecify  LEMIT  in  NAME1 
because  the  same  variable  is  set  equal  to  zero  in  NMEMSN  in  the  preceding 
part. 

The  energy  densities  obtained  from  this  example  are  compared  with  those 
of  Example  1  in  Table  9-3. 


9.  3  EXAMPLE  3 

We  now  consider  the  conical  cavity  illustrated  in  Figure  9-3.  Tnis  simple 
configuration  is  chosen  because  of  an  available  independent  solution.  The 
cavity  has  a  depth  of  10  inches  and  a  circular  opening  of  2  inches  radius. 

The  interior  surface  is  a  gray  body,  diffuse  emitting  and  reflecting,  with 
0.7  emiltance  and  at  a  uniform  temperature  of  300° K.  We  shall  calculate 
the  total  radiant  flux  out  of  the  cavity. 

The  input  data  for  this  calculation  is  listed  in  Table  9-4.  This  is  a  calcu¬ 
lation  of  internal  thermal  emission  rather  than  external  radiation  so  LEMIT 
is  zero,  as  given  by  default,  in  NAME1.  No  input  is  needed  in  NAME2 
because  neither  diffraction  nor  refraction  occurs  in  the  problem. 

The  conical  surface  is  divided  into  two  parts  for  importance  sampling.  The 
surface  definition  given  in  NMGEOM  is  clear  from  the  discussion  in 
Example  1.  No  special  constraint  needs  to  be  specified  in  NAME3. 
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TABLE  9-4.  INPUT  DATA  OF  EXAMPLE  3 


CONE  CAVITY  TO  VERIFY  EMISSION 
ENOTITLF 

1 

4NAMEI  NSEC=1 •  PMAP=2..  ISEXT=34 
4NAME2  4 

4NMGEOM  I  PLANE ( 1 >  =  3  *  ICONF <  1  >  =  1 *2* 
PTLCTR(l»l)=0»0»10tO«0*2*0»1.6« 

PTLCTR(l»2)=0f0»2»0«Q«0*l»fc»2» 

PTLCTR( l.jl =0*0.0* 0«0» 1«2» OS 
4NAME3  4 

4NAME<*  NSFACE  < 1 ) —  3  *  ISFACF U * 1 ) *1 »2* 1* 

KCOEFUf  1)«0«0«1«  ISFTP  ( 1 « 1 ) =5»5«2*  IMPSF  ( l  ,  1 )  *  1 . 0«  0  * 

NIMPS ( 1 )  =  1 «  IMPORT ( 1 « 1 ) =3  «  NH 1 NG ( 1 1 1 )  =  i «  N5EG<1*1)*5S 

SNMCOAT  C0AT(l.l)s.7»0*.3*  COAT ( 1 .2)=. 7.0* .31 
4NMEMSN  NP0LL=2«  LEMIT*200,  POLE  ( 1  ♦  1)  =0  *  0  *<♦  ♦  POLE U *2) *1  *  0 • 1 . 
NSFEM ( 1 ) si 1 1 1  KSECEM ( 1 « 1  )  =  1 »  KSECEM ( 1 »2) -1 « 

KS I  DEM (1*1)— 1*  KS IDEM ( 1 <  2 ) s2 «  TEMP ( 1 ) *300 * 300 . 

DELTAEs.001*  OEPRT=.l,  LPFN=8$ 
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FIGURE  9-3.  CONICAL  CAVITY 


The  values  given  to  NSFACE,  ISFACE,  KCOEF  and  ISFTP  are  clear  in 
reference  to  the  figure.  Radiant  energy  is  to  be  accumulated  at  the 
opening,  so  surface  3  is  made  the  important  surface.  While  the  rays  from 
surface  1  are  split  toward  the  important  surface,  those  from  surface  2 
are  allowed  to  select  random  directions.  The  latter  are  not  split  because 
they  are  so  close  to  the  important  surface  that  the  large  variation  of  the 
important  ray  energy  can  cause  significant  error;  see  Subsection  4.4. 
Furthermore,  since  the  rays  have  a  high  probability  of  selecting  the 
directions  toward  the  opening,  there  is  no  need  to  importance  sample. 

The  area  of  surface  3  is  comparable  to  the  square  of  the  shortest  distant 
from  surface  1  to  surface  3  so  it  is  advisable  to  partition  the  important 
surface.  It  is  divided  into  5  segments. 

Surfaces  1  and  2  have  an  absorption  of  0.7  and  a  diffuse  reflectance  of  0.3. 

Parameters  for  internal  emission  are  specified  in  NMEMSN.  Since  rays 
from  surface  2  have  a  high  probability  of  striking  the  exit,  more  rays 
should  be  generated  there  than  called  for  by  the  area  ratio.  This  is 
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achieved  by  placing  two  poles  in  the  system.  The  first  pole  is  placed  near 
the  center  (0,  0  4)  and  the  second  near  surface  2  (1,  0,  1).  The  latter  is 
placed  one  inch  off-axis  for  the  purpose  of  generating  rays  from  surface  2 
with  various  levels  of  energy. 

There  is  one  emitting  surface  associated  with  each  pole.  Both  surfaces  are 
in  section  1.  Surface  1  is  associated  with  pole  1  while  surface  2  with  pole 
2.  Both  interior  surfaces  are  at  the  uniform  temperature  of  300° K. 

The  distance  of  a  typical  pilot  ray,  the  ray  from  a  pole  to  a  corresponding 
emitting  point,  is  in  the  order  of  3  inches.  Assuming  0.7  for  a  typical 
value  of  cos  <f>,  see  Figure  3-2,  the  energy  of  a  ray  at  emission  is  in  the 
order  of 

2  4 

2n  s  e g  T 

(cos  <f|) 


or 


e 


2n  32  x  0,7  x  3.657  x  10~U  x  30Q4 
0.  7 


watts 


=  16.  7  watts 

_3 

The  cut-off  energy  level  of  10  watts  is  sufficiently  small  compared  with 
the  above  value.  The  ray  tracing  is  limited  to  8  reflections. 

Tracing  200  pilot  rays  from  each  pole,  we  obtain  an  emissive  pe'cr 
of  3.24  watts  at  the  exit,  giving  an  equivalent  cavity  emittance  of  0.871. 
This  compares  favorably  with  a  value  of  0,879  obtained  from  an  independent 
source  (Reference  3). 


!).  4  EXAMPLE  4 

The  conical  cavity  in  Example  3  was  chosen  so  that  comparison  can  be 
made  with  the  result  of  an  independent  solution.  We  will  now  proceed  to 
analyze  a  cavity  and  demonstrate  that  the  GUERAP  program  can  be  used 
in  solving  a  wide  variety  of  problems  on  thermal  radiation  of  which  solutions 
are  heretofore  unavailable. 

As  shown  in  Figure  9-4,  the  cavity  is  a  circular  cone  of  10-inch  depth  and 
22.5  degrees  half-angle.  The  front  end  is  covered  by  a  plane  with  a  perfect 
circular  aperture  of  2,5  inches  radius.  Both  interior  surfaces  are  gray 
bodies  with  0.7  emittance  and  at  a  uniform  temperature  of  300°  K.  As  in 
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FIGURE  9-4.  MODIFIED  CONIC  AL  CAVITY 


the  last  example,  we  are  looking  for  the  total  raaianl  energy  leaving  the 
cavity  and  the  equivalent  emittancc  at  the  exist  aperture. 

The  input  data  given  in  Table  9-5  are  rather  straightforward.  The  exit 
aperture  is  the  important  surface  for  surface  1  while  rays  from  surface  2 
are  allowed  to  be  randomly  distributed. 

No  particular  emphasis  seems  necessary  for  emission  from  any  part  of 
the  interior  and  a  pole  in  the  middle  of  the  system  is  sufficient.  There  are 
two  surfaces,  surfaces  1  and  2,  associated  with  the  pole  and  both  are  in 
section  1 . 

The  cut-off  energy  level  is  set  at  10'^,  substantially  lower  than  necessary. 
However,  this  will  not  significantly  prolong  the  calculations  because,  all 
the  rays  are  limited  to  eight  reflections. 

With  100  pilot  rays,  the  result  shows  that  energy  is  leaving  the  exit 
aperture  at  a  rate  of  4.58  watts,  giving  an  equivalent  opening  emittancc  of 
0.788.  The  latter  is  substantially  lower  than  that  of  the  last  example 
because  cf  higher  ratio  of  the  opening  radius  to  the  cavity  depth. 
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TABLE  9-5.  INPUT  DATA  OF  EXAMPLE  4 


CONE  CAVITY  WITH  REDUCED  EXIT  OPENING 
END  TITLE 

1 

JNAMEl  NSEC= 1 ♦  RM AP=2 . 5 ♦  ISEXTOS 
SNAME2  * 

SnMGEOM  IPLANE  ( l  ICOmE<1>=1. 

PTlCTPI 1,1 >*0*0«10«0* 0.0 *0*4. 1421. 

PTLCT«( 1,2> =0.0. 0.0*0* 1*4.1421.2.5* 

PTLCTPI 1.3) =0.0. 0.0.0. 1*2.5.04 
tNAMEJ  i 

SNAME4  NSFACE(I)=3.  I  SPACE  (l*l)sl«2*3«  KCOEF ( 1 • 1) *0  *  1  •  1 » 
ISFTP(1«1>*5«5*2*  IMPSP ( 1 » 1 ) =1 .0.0.  NlMPS ( 1 ) -1  * 
IMPOPT(1,1)s3«  NR IMG (  1  .  1 )  =  l  ♦  NSF.G(l,l)sli 
SNMCOAT  COAT ( 1 « ) )=.7«0«.3.  COAT ( 1 . 2) = . 7 » 0 ♦ . 31 
tNMEMSN  NPOLEsl*  LEMITslOO,  POLE  (1*1)  *0*0***.  NSFEM(l>*2. 
KSECEM  ( 1  *  1 ) S1 « 1 «  KSIDEM1.1  )*1*2*  TEMPI  1 ) =300.300. 
UELTAE=l.t-7.  OEPRT=l.E-7*  LRFN=B$ 
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9.  5  EXAMPLE  5 

We  now  consider  an  optical  system  illustrated  in  Figure  9-5,  The  input 
data  is  simplified  by  the  fact  that  the  system  is  axially  symmetric  around 
the  Z-axis.  Sun-light  with  a  quarter  degree  divergence  enters  the  system 
at  45  degrees  from  the  axis.  The  first  six  sections  consist  of  a  series  of 
baffles  aimed  at  attenuating  the  off-axis  energy.  From  the  seventh  section 
on  there  are  an  aspheric  corrector  lens,  two  meniscus  lenses  and  one 
spherical  mirror  on  the  right  hand  end.  The  on-axis  energy  is  specularly 
reflected  by  the  mirro”  and  focused  onto  the  detector,  surface  56. 

All  the  required  dimensions  are  clearly  indicated  in  the  figure.  Also  shown 
are  surface  numbers,  numbers  in  small  circles.  Most  of  them  are  either 
real  surfaces  or  apertures,  while  some  are  imaginary  surfaces  used  ex¬ 
clusively  for  importance  sampling.  We  will  first  briefly  describe  the  sur¬ 
faces  in  each  of  the  sections.  It  should  be  emphasized  here  that  boundaries 
of  sections  are  chosen  chiefly  for  convenience  in  the  handling  of  data  and 
efficiency  in  the  execution. 

Section  1  is  bounded  by  a  cylinder,  a  plane  baffle,  a  conical  baffle,  and 
two  apertures.  Referring  to  Figures  2-4  and  2-5,  we  see  that  if  we  assign 
7  and  11  to  the  surface  numbers  of  the  baffle  knife  edges,  the  baffles  will 
consist  of  surfaces  7  through  14.  All  of  these  surfaces  except  9,  12  and 
14  are  invluded  in  the  section.  Surfaces  49  and  50  are,  respectively,  the 
entrance  and  the  exit  apertures,  of  the  section.  In  addition  we  place  a  plane 
annular  ring  with  a  larger  inner  constraint,  surface  46,  directly  in  front 
of  surface  10.  The  purpose  of  this  surface,  which  is  to  be  used  for  impor¬ 
tance  sampling,  will  become  clear  later. 

Sections  2,  3,  5  and  6  are  essentially  similar  to  section  1  except  that  no 
additional  surface  like  46  is  created. 

Section  4  is  bounded  by  a  cone,  two  conical  baffles  and  two  apertures. 
Similar  to  section  1,  we  place  a  cone  with  larger  inner  constraint,  surface 
47,  directly  in  front  of  surface  22. 

The  small  circular  cylinder  in  the  middle  of  the  system  is  black  coated 
on  all  sides.  Surfaces  of  section  7  thus  include  cylinders  4,  41,  and  43, 
planes  35,  42  and  43,  baffle  surfaces  31,  32  and  34,  and  aperture  55. 

Section  8  is  the  interior  of  the  aspheric  corrector  bounded  by  planes  35 
and  36,  and  cylinders  41  and  43,  Section  9  is  the  space  between  the  first 
two  lenses.  It  is  bounded  by  planes  36  and  42,  cylinders  41  and  43,  and 
sphere  37. 
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Section  10  is  the  interior  of  the  first  meniscus  lens  bounded  by  plane  42, 
cylinders  5  and  43,  and  spheres  37  and  38.  Section  11  is  the  interior  of 
the  second  meniscus  lens  bounded  by  cylinders  5  and  43,  and  spheres  38 
and  39. 

Section  12  is  the  space  between  the  second  meniscus  lens  and  the  spherical 
mirror.  It  is  bounded  by  planes  45  and  56,  cylinders  5  and  48,  cone  6, 
and  spheres  39  and  40.  Surface  56  is  a  circular  detector,  where  energy  is 
to  be  collected,  placed  inside  the  annular  ring  45.  Surfaces  5  and  48  are 
the  same  cylinder  divided  by  a  plane  constraint  at  Z  =  4.5  inches  for 
importance  sampling  purposes. 

In  addition,  there  are  seven  more  planes  to  be  used  as  important  surfaces. 
57  is  on  the  first  baffle.  In  the  vicinity  of  the  fourth  baffle,  there  are 
planes  59,  60  and  61.  In  section  7,  there  is  the  annular  ring  58  surrounding 
circle  44.  In  section  12,  we  have  62  approximately  coinciding  with  the 
mirror  and  63  right  on  the  detector.  Dimensions  for  these  surfaces  are 
not  clearly  indicated  in  the  figure.  However,  they  can  be  easily  resolved 
from  the  input  data,  shown  in  Table  9-6. 

Sun-light  with  a  quarter  degree  divergence  comes  into  the  system  at  45 
degrees  from  the  normal.  So  we  have  SUNAGD  =  45  and  SUNDVD  =  0.25. 
Since  the  system  is  axially  symmetrical,  the  incident  beam  can  be  in  any 
direction  and  the  azimuthal  angle  SUNAZD  assumes  the  default  value  zero. 

To  estimate  the  minimum  energy  level  DELTAE,  we  first  assume  an 
extremely  small  value  and  generate  a  small  number  of  rays  with  reason¬ 
able  number  of  interceptions  LRFN.  The  result  of  the  above  trial  compu¬ 
tation  will  show  the  desirable  level  of  DELTAE. 

Surfaces  49  and  56  are,  respectively,  the  entrance  and  the  exit  apertures. 

Plane  surfaces  are  35,  36,  42,  44,  45,  46  and  49  through  63.  The  inside 
radius  of  42  is  obtained  from  the  intersection  of  the  surface  with  sphere  37. 
Enough  accuracy  should  be  maintained  in  this  value  to  assure  that  no  gap 
exists  between  the  two  surfaces.  It  is  advisable  to  use  surface  37  as  the 
inner  constraint  of  surface  42. 

Plane  46  of  inner  constraint  2. 3  inches  is  placed  closely  in  front  of  plane 
10.  It  will  thus  intercept  all  the  rays  that  come  within  its  boundary.  This 
will  allow  the  command  to  follow  the  instruction  concerning  importance 
sampling  under  surface  46  instead  of  10.  Surfaces  49  through  55  are  a 
series  of  circular  apertures.  Their  inside  radii  are  all  zero. 


9-17 


:t:u»(>-ooH 


TABI.K  9-B.  INPUT  DATA  OK  KX  AMPLE  5 


EAAMRLf  o 
ENDTITLE 

1 

ShAME  l  5t)NA6D  =  *b..bUN0vns.;>b.l)tLUf.  =  l  .t-15*DtPR!*i.E-l2*LCM!!*bO. 
LKFN=.h)»  lSEtL>R=10.NPRlNT  =  3nO,NERR=100.NbECM2*R4AP«.  19375.KPR|NT*0* 
ISEnT=49, IStAT=b65 

VJAME2  RINOEa<  nsl«  I  *1.1. 1,1.1  .l.4o 18, 1.1.025.1,625.  IS 

SNMC.EOM  IPLA  JF  (  1  )=.3b«3o«42. 44.45.44. 49. 50, :>1*52»S3»S4.55*56.S7.56. 

b9»f>0«61 *  62  »  6 1 » 

PTLCTR!  1  ,35)=2"0.  .6,655,2*0.,  1  ..  1.5936 ..45* 

PTLCTW(  I  ,36)  *?•<).  '<5.155.2*0.  .  1  •«  1  .5936. .45* 

PTLCTR ( 1 .42)  ■-•2*0., 4. Obi  .2*0.,  I  .*2. 002*.  I  .54854, 
PTLCTR(l.44)=2*0..7.906»2*0.»l.»,4b,0.» 

PTLCT4(  I  •  4S )  =  2°0  •  »4 .624*.  ,2°0.*1.,  .45,  ,19375. 

PTLCT* ( 1 .4b) =2*0. 1 1 7.2*099.2*0., 1 . . 4 . a 7b . 2 . 3 • 
PTLCTR(l,49)=2*Q.,17.33»2*0.,l.»l .645,0  •  , 
PTLCTR<1.50)=2*G..16.46.2«0..1.,J.005,0.» 
PTLCTR(l»5l)=2*0.»15.36,2*0.»l.,2.315.0.. 

PTLCTR (l.b2) =2*0..  14. 14R23, 2*0.. i..l  .05.0.. 

PTLCTR (1,53) =2*0..  11  .46472.2*0. »1..> .44.0.* 

PTLCTR  ( 1  ,54)=2*U.,9.58»2*Q,.1,»1»375»0.* 

PTLCTR (1,55) =2*0. *9. 03. 2*0.,  l.,l. 375.0. • 

PTLCTR( 1, 56) =2*0.. 4. 6294, 2*0.. 1.*. 19375.0.* 
pTLCT9( 1.57) *0.0.1 7.33*0. 0,-1, 2. 3. 1.645* 

PTLCTR ( 1.S8) *0.0 *7. 905. 0.0, 1.1. 4*. 45* 

PTLCT9 11.59) *0.0. 15. 38.0.0, 1,2. 3 15.0. 

PTLCTR ( 1 .60) *0,0. 14.14923.0,0, 1,1 .535.0. 

PTLCTR( 1.61) =0.0.1 3.9. 0.0, -1.2.47 .1.535. 

PTLCTR(1. 62 )*0» 0,0, 0,0,1, 1,76.0. 

PTLCTR (1.63 )*C. 0,4. 6294,0,0,-1.. 19375.0* 

I CONE (1)*1, 2, 3*4*5.6, 4J.43.47. 4H* 

PTLCTR ( 1,1) *0.0. 17.33*0.0* 13. 5* 2*4. 87*. 

PTLCTR ( 1 .2) *0.0. 13.5,0,0. 1 1 .04,4.875*3.625* 

PTLCTR (1,3) *0.0. II.  04. 0.0. 9. 58. 2*3. 625* 

PTLCTR (1.4) *0.0 *9. 58* 0.0 *6. 061 *2*2. 0025* 

PTLCTR (1.5 ) *0.0 *6. 06 1*0. 0*4. 5* 2*1. 77. 

PTLCTR (1.6) *0.0*3. 118.0*0*. 2418. 1.77*2.0 72* 

PTLCTR ( 1,41) *0.0. 6. 655*0 .0,6. 06 1*2* 1.5936. 

PTLCTR ( 1.43) *0.0* 7. 905, 0,0, 4. 6294. 2*. 45* 

PTLCTR ( 1.47) *0*0*1 3.9 1657, C *0 . 1 3.49990*2.4 7*4,875* 

PTLCTR ( 1,48) *0.0. 4. 5*0. 0*3. 11 8, 2*1.7 7* 

ISPmER(1)*37. 38. 39*40. 

PTLCTR (1.37) *0,0, 9, 0,0.5.678*1.54864,. 45. 

PTLCTR (1,38) *0*0, 9. 638. 0.0, 5. 278, 1.77,. 45* 

PTLCTR( 1 .39) *0.0, 9, 0.0. 4. 6063* 1.77*. 4R. 

PTLCTR ( 1,40) *0.0, 9. 0.0* 0,2. 072.0, 
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TABLE  9-6.  INPUT  DATA  OF  EXAMPLE  5 

IBAFF (1 ) =7  *  1 1  *  15 , 19*23*27, 3i * 

PTLCT9 (1  *  7 ) *0  *  0 , 1 7 • 33* 0 .0. I , 1 .645, 0 , 

PTLCTP( 1* 11 ) *0.0 *lh. 46.0. 0.1 *3. 005.11. 
PTLCTW(1*15)=0*0.15.38*0*0, 1.2.315*  11, 

PTLCT9(1*19)=0* 0*14.1 4923*0 *0*1* 1.535* 11, 

PTLCTU (1.23) =0*0. 11. 46472. 0,0. 1,1. 44*1 1* 

PTLCT9( 1.27) *0*0.9. 58. 0,0. 1*1.375,0. 

PTLCTWU,  31)  *0*0. 9. 03*0  *0*1  *1.375*0* 
RKNIFE=.00003*TBAFF=.05*DBAFF=.001 ,ABaFFD=15* 

4NAME3  NST9T  d)=0*0*0*0*2*l  *2*2*0*1*2«2*1*1*2*2*1.1*2*2* 
NSTPT(21)  =  l.l*2*2*l*l*2*2*l*l»2*2*l,l*2*0*2*Q*l)*Q* 

NSTPT (4l)*l* 2*2*14 

tNAME4  NSFACE (1) =9 *8, 8*9*8,8*10*4*5*5, 4*7* 

ISF ACE (1* l ) *1*7 *8, 10, 11* 13. 46*49*50. 

I  SPACE ( 1.2) *1*1 1*12* 14, 15. 17. 50*51* 

ISF ACE (1,3) *1*15* 16* IB, 19.21 *51 *52. 

ISFACE (1.4) =2. 19*20, 22,23,25. 47.52.53, 

1SF ACE ( 1  *  5 ) *3  *  23  *  24  *  26*27  *  29*  53*54* 

ISFACE (1,6) =4. 27, 28, 30, 3 1*3 3* 54, 55, 

I  SPACE (1*7) *4,31  *32*34*35*41*42*43*44.55* 

ISFACE (1*8) =35*36*41,43* 

ISFACE (1*9) *36* 37, 41,42*43, 

I  SPACE (1,10) =37* 36*42*5*43, 

ISFACE (1*11)* 38 *39, 5, 43* 

I  SPACE (1,12)* 39, 40*5. 6, 45 *48* 56* 

KCOEF (1*1)=0*0*0*1*0*1. 0*0.1, 

KCOEF (1,2)=0,0*U*0*0*1*0,1* 

KCOEF (1.3) =0.0.0 *0*0, 1,0,1, 

KCOEF  ( 1.4) =0.0*0*0*0* 1,0*0*  1, 

KCOEF (1,5) *0.0, 0*0*0* 1,0,1, 

KCOEF (1*6)=0*0,0*1*0*1*0«1* 

KCOEF (i*7)=0*0*0,l*l*l*l*l,l,0* 

KCOEF (1,8)*0«1,0,1, 

KCOEFO  *9)*0*0, 0*1*1* 

KCOEF (1.10)*1,0, 0.0.1* 

KCOEFd.l  1  )*1, 0*0,1. 

KCOEF (1. 12) =1,0.0. 0*0*0. 0, 

ISFTP(1,1)=7*5*1,2* 

ISFTP(l,2)*b*5*l»2* 

ISFTP(1*3)=6*5,1.2, 

ISFTP(l*4)*7*5*1.2* 

ISFTPd  ,5)*6*5*1,2* 

ISFTP(1,6)*6*5*1*2* 

ISFTPd.  7)  =4*5. 4, 4*5*1, 

ISFTP(1*8)=3.4*5.5* 

ISFTP(1*9)=3»4*5*5*5* 
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TABLE  J-6.  INPUT  DATA  OF  EXAMPLE  5 

ISFTP(1*10)*3*4«5*5*5* 

ISFTP(1*11)*3****5*S* 

ISFTP(I*12)s3«S*b»?* 

lM°SF(l*l)sl*2»2*2*l*l*0*2.1* 

1MPSF<1*2)*0«1*1*1*2*2*1«0, 

IMPSF (1*3) *0*1, 1*1*0* 0.1.0, 

IMPSF (l«4)3l*2*2*2«l*l*0*2,l* 

IMPSF ( 1  *5)  s0 *  1  • 1 • 1 *0* 0 « 1 *0. 

IMPSF ( 1 *6)s0* 1 * 1* 1.0. 0.1.0. 

I MPSF (1«7)*0*1*1*1«0*0*0*0*0*1* 

IMPSF  (1*12)*1«2*1*4*1*4*0* 

NIMPS(l)s2»2»l*2»l*l*l*0*0*0*0*2* 

IMPORT  ( 1 « 1 )  =S7«Sh*  NpINGU*  1  )*1.3*  NSEG  (1  ♦  1 )  *4*  1 * 

IMPORT ( 1 *2) =69*57*  N9lNG( 1 *?) =1 . 1 «  NSEG ( 1 •?) *4*4 » 

IMPORT  ( I «  3)  =60  *  NP  ING  ( 1  «  3)  =  1 «  NSEG(l*3)*4. 

lMP0RT(l«4)=6l*b8»  NR ING ( 1 *4 ) = 1 , 3*  N$EG(  I  *4)  *4. 1  « 

IMPORT  (1*5)  =54*  NPINGU  *S)*1.  NSEGU«b)*4* 

IMPORT  ( I  *6)  *55*  NRlNG(l*6)*l.  NSEG(1*6)*4* 

IMPORT  <  1  ,7)=35*  NRlNGU  *7)=1.  NSEG(1«7)*4* 

IMPORT (1*12)=62*63«  NRING( 1  *  12) =1 « l  •  NSEGl 1 • 12) *4* IS 

SNMCOAT  COAT ( 1  *  1 )  =  .96*0* .04,  COAT <1 *2 > =.96 » 0 ♦ .04*  COAT  1 1  *3)  *.96*0*  .04* 
COAT(1»4)=.96»0».04»  COAT ( 1  * S )  =  .96* 0 . . 04 •  COAT ( 1 *6) *.96*0* .04* 

COAT (1*  7)  =  .96«0*.04»  COAT ( 1  * M )  =  .96* 0 • .04*  COAT (1 *9)*. 96*0* *04* 

COAT  ( 1 « 10)  =  .96*0*  .04,  COAT  ( 1  *  1 1 )  =.96.0*  .0*»*  COAT  (1*12>*.96*0*.04. 

COAT  ( 1 • 1 3 ) =  *96  «  0 • «  04 •  COAT  (l«l4)=.9to.0*,04,  COAT  (1*15)*. 96 *0*. 04* 

COAT  ( 1 . 16)  =  .96*0*  .04*  coat  (1*17)*.  96*0*.  0***  COAT (1 • 1H) *.96* 0* .04 • 

COAT ( 1 « 19)=.96«0* .04,  COAT ( 1 *20) a. 96.0* .0* *  COAT ( 1 *21 ) *.96*0* .04* 

COAT (1*2?)=. 96*0 *.04,  COATU«23)s.96.0*.0**  COAT  (1  *24)  *.96*0*  .04* 

COAT  ( l  ,25)  =  .96*0*  .04*  COAT  ( 1  *26)  *  .96, 0 *  .04 •  COAT  ( 1  *27)  *.96*0 •  .04 • 

COAT  (l,2fl)=.96,0».04,  COAT < 1 *29) =.96,0* .04,  COAT  ( 1  *30 ) *.96*0*  .04 . 

COAT  U, 31)  =  .  96*0*.  04*  COAT  ( 1  *32)  =  .96.0*  .0**.  COAT  U  *33)  *.96*0*  .04* 

COAT (1,34) =.96*0 *.04,  COAT  ( 1 « 4 1 )  =  . 96. 0 *  •  0***  COAT  ( 1  *42) *,96*0* .04, 

COAT  (1*43)3. 96  *0*.04«  COAT ( 1 *44 ) = , 96, 0 • .04,  COAT  ( 1  *45)  *.96*0  »  .04* 

COAT  ( 1 ,46)  s.96*0*  .04*  COAT  (1*47)*, 96,0*. 04,  COAT  U  *40)  *.96*0*  .04* 

COAT (1.3S) =. Cl*. 02*. 01 *0*0.. 95*. 01*  COAT (1  * 36) *.01  * .02* .0 1 *0*0 • .95* .0 1 • 
COAT (1,37)  =  .01* .02* .01 *0*0.. 95* .01*  COAT ( l • 36) *.01  * .02* .01 *0*0* .95* .01 • 
COAT ( 1 . 39) =.01 • .02* .01 *0*0* .9S* .01 ♦  COAT ( 1 *40 ) *.0 1 • .97* .0 1 • .01  * 
CQAT(14*40)sb9.S9* 

SNMFMSN  LF.M  I  T  =  0S 
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Cone  and  cylinder  surfaces  are  1  through  6,  41,  43,  47  and  48.  Sphere 
surfaces  are  37  through  40.  The  parameters  of  the  above  surfaces  are 
obvious  with  reference  to  Figures  2-1  through  2-6. 

All  the  baffles  have  an  edge  radius  of  0.00003  inch,  a  thickness  of  0.05 
inch,  and  a  blade  angle  of  15  degrees.  For  the  plane  baffles,  we  assume 
a  halfwidth  of  0.001  inch  for  the  hyperboloid;  see  Figure  2-5. 

At  mentioned  in  Section  2,  there  are  four  surfaces  for  each  baffle.  Param¬ 
eters  are  read  in  for  the  knife  edge  surfaces  7,  11,  15,  19,  23,  27  and  31. 

7  27,  and  31  are  plane  baffles  and  the  rest  are  conical.  These  are 
indicated  by  the  last  argument  of  each  group.  As  is  given,  the  inclination 
angles  for  the  conical  baffles  are  11  degrees. 

The  standard  constraints  given  by  the  calculation  of  baBic  configurations 
are  satisfactory  for  most  surfaces.  However,  slight  modifications  of 
certain  constraints  are  sometimes  necessary.  Each  of  the  conical  sur¬ 
faces  1  through  4  has  two  standard  constraints.  It  is  clear,  however,  that 
these  cones  can  be  allowed  to  extend  indefinitely  and  the  rays  in  a  section 
will  not  be  falsely  intercepted  by  the  extended  portion  of  the  cones.  The 
rays  will  be  intercepted  by  other  surfaces  because  these  are  at  closer 
distances.  Thus,  we  can  remove  all  the  constraints  from  these  surfaces. 
This  is  achieved  by  specifying  NSTRT(I)  =  0  for  I  *  1,  2,  3  and  4. 

Similarly,  no  constraint  is  required  on  the  right  hand  side  of  conical  sur¬ 
face  6.  A  constraint  is  required  for  the  left  hand  side,  however,  because 
the  conical  surface  would  have  extended  into  the  section.  Referring  to 
Paragraph  2.1.2  and  the  array  PTLCTR  for  surface  6,  we  see  that  the 
fl  'at  standard  constraint  is  a  plane  on  the  left  at  Z  >  3.118  while  the 
second  one  is  a  plane  on  the  right  at  Z  *  0.2418.  To  eliminate  the  second 
constraint  while  keeping  the  first  one,  it  is  necessary  only  to  specify 
NSTRT(6)  *  1. 

Referring  to  Subsection  2.1.5,  we  see  that  surface  10,  surface  4  of 
Figure  2-5,  is  constrained  only  by  the  knife  edge.  It  is  clear  that  the 
surface  can  be  allowed  to  extend  indefinitely  on  the  outside.  Thus  the 
number  of  constraints  remains  at  1. 

Surface  numbers  for  each  section  are  clearly  indicated  in  Figure  9-5. 
Considering  section  7.  we  have  10  surfaces,  that  is  4,  31,  32,  34,  35,  41, 
42,  43,  44  and  55. 
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The  arrays  KCOEF  and  ISFTP  are  arranged  in  the  same  way  as  the  surface 
numbers.  We  will  again  consider  section  7.  Surface  4  is  a  cylinder  with 
the  region  inside  so  KCOEF  is  zero.  Surface  31  is  a  hyperboloid  with  the 
region  on  the  side  including  the  center  so  KCOEF  is  zero.  (See  Figure  2-7.) 

Surface  35  is  a  transmitting  lens  toward  the  exit  so  ISFTP  is  4.  Surface 
55  is  an  aperture  toward  the  entrance  so  ISFTP  is  1.  The  rest  of  the 
surfaces  in  the  section  are  all  non -transmitting  surfaces  so  ISFTP  is  5  for 
each  of  them. 

Important  surfaces  are  5^  and  58  for  section  1,  59  and  57  for  section  2, 
etc.  The  first  important  surface  of  section  1,  surface  57,  is  divided  into 
five  rings  and  three  segments,  etc. 

An  understanding  of  the  characteristics  of  the  system  will  assist  in  the 
selection  of  important  surfaces.  The  incident  beam  45  degrees  off-axis  is 
intercepted  by  surfaces  15,  17  or  21.  We  are  concerned  with  the  energy 
that  after  many  transmissions  and  reflections  eventually  reaches  the 
detector.  As  it  appears  then,  we  will  be  interested  in  the  energy  that  goes 
into  the  aperture  toward  the  detector.  For  example,  if  we  use  surface  50 
as  the  important  surface  for  surfaces  7,  8,  10  and  49,  we  will  obtain  the 
energy  entering  section  2.  By  repeating  the  process  for  all  the  sections, 
we  will  arrive  at  the  energy  that  reaches  the  detector. 

As  we  learned  from  experience,  however,  we  do  not  obtain  the  best  result 
using  this  procedure.  For  example,  when  a  ray  is  reflected  off  surface  10, 
the  portion  of  the  energy  going  through  aperture  50  is  rati.er  small.  After 
so  many  reflections  and  splitting  the  energy  for  the  ray  is  so  small  when 
it  reaches  the  detector  that  it  is  insignificant  when  compared  with  the  total 
energy  at  the  detector.  Thus,  automatic  choice  of  exit  aperture  of  a  section 
as  the  important  surface  is  not  necessarily  the  best  approach. 

Consider  a  ray  being  reflected  off  surface  17.  The  "shortest"  way  for 
this  ray  to  reach  the  detector  is  by  a  single  reflection  off  surface  10  and 
straight  down  through  the  lenses.  The  energy  of  this  ray  will  be  reduced 
by  importance  sampling  splitting  at  surfaces  17  and  10,  and  absorption 
and  scattering  at  the  lenses  and  the  mirror.  Because  the  ray  will  hit  the 
mirror  oif-axis,  the  energy  that  reaches  the  detector  will  be  further 
reduced  by  the  mirror  scattering  coefficient.  However,  the  resulting 
energy  of  such  a  ray  can  be  significantly  higher  than  the  one  described 
above. 

Considering  section  1  then,  we  chose  the  first  Important  surface,  surface 
57,  as  the  important  surface  for  surfaces  1,  11,  13  and  50.  Thus,  every 
ray  from  any  of  these  surfaces  will  be  split  toward  surfac  *  10. 
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Thr  second  important  surface  of  section  1  is  the  ring  58,  in  front  of  the 
aspheric  collector.  By  assigning  this  as  the  important  surface  for  surfaces 
7,  8,  10  and  49,  we  focus  our  attention  on  the  energy  that  is  going  toward 
the  lenses.  Note  that  surface  10  occupies  only  the  middle  portion  of  the 
baffle  surface.  When  a  ray  hits  the  baffle  in  the  outer  area,  Uie  surface 
number  will  be  46  and  there  will  be  no  splitting.  Hence,  the  values  of 
IMPSF  are  0  for  surface  46,  1  for  surfaces  1,  11,  13  and  50,  and  2  for 
surfaces  7,  8,  10  and  49.  , 

There  are  three  kinds  of  surfaces  in  terms  of  radiation  properties. 

Surfaces  1  through  34  and  41  through  48  are  black  coated  structural 
surfaces.  They  have  a  96  percent  absorptance  and  a  4  percent  diffuse 
reflectance.  Surfaces  35  through  39,  the  lens  surfaces,  have  1  percent 
absorptance,  ?  percent  specular  reflectance,  1  percent  diffuse  reflectance, 
95  percent  specular  transmittance  and  1  percent  diffuse  forward  scattering 
coefficient.  The  mirror,  surface  40,  has  1  percent  absorptance,  97  per¬ 
cent  specular  reflectance,  1  percent  diffuse  reflectance,  and  1  percent 
near  specular  reflectance.  The  angle  constant  of  the  exponential  decay  for 
the  near  specular  reflectance  is  obtained  from  an  analysis  of  the  mirror 
scattering  measurement. 

The  user  is  warned  that  an  execution  of  this  example  requires  a  relatively 
large  amount  of  computing  time  due  to  its  complexity:  approximately  15 
minutes  on  a  CDC-6600  computer.  This  problem  was  taken  from  an 
actual  analysis  performed  during  the  development  stage  of  the  GUERAP 
program.  It  is  included  to  show  the  procedure  for  setting  up  an  input  data 
set  for  a  practical  problem. 


9-23 


3366-008 


REFERENCES 


1.  "Use  of  Expected  Value  Techniques  to  Reduce  Computer  Time 
Requirements  for  Multi-stage  Baffle /Optical  Systems  '  Appendix  1, 
Unwanted  Energy  Rejection  Technique  Development  (UERTD)  Progress 
Report  No.  2,  10  September  1971,  Honeywell  Document  3366-PR2. 

2.  ibid,  "Diffraction  Modeling  ",  Appendix  2. 

3.  "GUERAP  User's  Manual",  T.  S.  Chou  and  B.K.  Likeness,  6  April 
1972,  Honeywell  Document  UERTD  TM  3366-007. 

4.  "Fundamentals  of  Optics",  Jenkins  and  White,  McGraw-Hill,  1957. 

5.  "Principles  of  Optics  ",  Born  and  Wolf,  Pergamon  Press,  1959. 


9-24 


